Carrega de les llibreries que es necessiten

library(arules)
library(C50)
library(caret)
library(class)
library(cluster)
library(dplyr)
library(factoextra)
library(FactoMineR)
library(fpc)
library(ggpubr)
library(ggplot2)
library(grid)
library(gridExtra)
library(mclust)
library(plyr)
library(purrr)
library(randomForest)
library(rattle)
library(RColorBrewer)
library(rpart.plot)
library(rpart)
library(vcd)

1 Introducció


1.1 Presentació

Aquesta pràctica cobreix de forma transversal l’assignatura.

Les Pràctiques 1 i 2 de l’assignatura es plantegen d’una forma conjunta de manera que la Pràctica 2 serà continuació de la 1.

L’objectiu global de les dues pràctiques consisteix en seleccionar un o diversos jocs de dades, realitzar les tasques de preparació i anàlisi exploratòria amb l’objectiu de disposar de dades llestos per aplicar algoritmes de clustering, associació i classificació.

1.2 Datasets

A partir de la pràctica 1 vam generar des de 2 datasets originals (retail i bank) 5 datasets resultants que farems servir per als diferents exercicis de la pràctica 2.

  • bank-clean-discret: conjunt de dades obtingut a partir del dataset de marqueting bancari on s’han estat discretitzats alguns camps, eliminat alguns camps nulls i eliminar outliers.

  • bank-discret: conjunt de dades obtingut a partir del dataset de marqueting bancari on s’han estat discretitzats alguns camps. No s’han eliminat nulls ni outliers.

  • retail-customer: conjunt de dades obtingut del dataset de venta online on s’han generat nous atribtuts per a tenir dades de comportament dels clients.

  • retail-product: conjunt de dades obtingut del dataset de venta online on s’han generat nous atribtuts per a tenir dades dels productes.

  • retail-association: conjunt de dades obtingut del dataset de venta online on s’han generat transaccions de les compras realitzades.


2 Associació: Generació de regles a partir de Regles d’associació.


A partir del dataset que hem preparat en la pràctica 1 anem a generar regles d’associació.

2.1 Carreguem el fitxer

retail.association <- read.csv('retail-association.csv', sep = ',')
summary(retail.association)
##    InvoiceNo        StockCode                                  Description    
##  Min.   :536365   85123A :  1686   WHITE HANGING HEART T-LIGHT HOLDER:  1680  
##  1st Qu.:549316   85099B :  1329   JUMBO BAG RED RETROSPOT           :  1329  
##  Median :562213   47566  :  1275   PARTY BUNTING                     :  1275  
##  Mean   :560833   20725  :  1207   LUNCH BAG RED RETROSPOT           :  1207  
##  3rd Qu.:572302   84879  :  1159   ASSORTED COLOUR BIRD ORNAMENT     :  1159  
##  Max.   :581587   22720  :  1120   SET OF 3 CAKE TINS PANTRY DESIGN  :  1120  
##                   (Other):330406   (Other)                           :330412  
##     Quantity                InvoiceDate       UnitPrice       CustomerID   
##  Min.   : 1.000   11/14/2011 15:27:   460   Min.   :0.000   Min.   :12347  
##  1st Qu.: 2.000   11/28/2011 15:54:   450   1st Qu.:1.250   1st Qu.:13994  
##  Median : 6.000   12/5/2011 17:17 :   442   Median :1.650   Median :15245  
##  Mean   : 7.477   10/31/2011 14:09:   384   Mean   :2.192   Mean   :15326  
##  3rd Qu.:12.000   11/23/2011 13:39:   370   3rd Qu.:2.950   3rd Qu.:16818  
##  Max.   :27.000   9/21/2011 14:40 :   370   Max.   :7.500   Max.   :18287  
##                   (Other)         :335706                                  
##            Country       Cancellation      TotalPrice    
##  United Kingdom:305151   Mode :logical   Min.   :  0.00  
##  Germany       :  7465   FALSE:338182    1st Qu.:  3.80  
##  France        :  6905                   Median : 10.08  
##  EIRE          :  5451                   Mean   : 12.77  
##  Spain         :  2046                   3rd Qu.: 17.40  
##  Belgium       :  1660                   Max.   :178.80  
##  (Other)       :  9504

2.2 Obtenció de transaccions

Fem una transformació per obtenir les transaccions de les compres segons el InvoiceNo. El camps que ens interessa veure és Description.

transations_split <- split(x=retail.association[,"Description"],f=retail.association$InvoiceNo)

# trasformem el data com a transactions. Cada transacció es una compra (InvioceNo) i els productes que ha comprat
transations <- as(transations_split,"transactions")
## Warning in asMethod(object): removing duplicated items in transactions
itemFrequencyPlot(transations, support = 0.04, cex.names = .6, col = rainbow(15), topN=30)

Mostrem els productes que més es repeteixen. En concret estem mostrant aquells que tenen un suport o freqüència >= 0.04 Això vold ir que, si tenim 16835 transaccions, estem mostrant x >= 16835 * 0.04, és a dir aquells que apareixen com a mínim 673,4 vegades.

Si llancem l’algoritme “apriori”, generarem directament un set de regles amb diferent support, confidence i lift.

  • El support indica quantes vegades s’han trovat les regles {lsh => rhs} en el dataset, com més alt millor. Es la “popularitat” d’un conjunt d’elements del dataset. On {lsh => rhs} indica que si es compra lsh, compra rhs.

  • La confidence ens parla de la probabilitat de comprar {rhs} si es compra {lhs} (lhs => rhs / lhs).

  • El lift és un paràmetre que ens indica quan d’aleatorietat hi ha a les regles. Un lift d’1 o menys ens indica que la regla és completament fruit de l’atzar. El lift ens diu eb una regla, com es produeix la regla en funció dels elements de la regla (support(lhs => rhs) / support(lhs) * support(rhs)).

Generem les regles per un support mínim de 0.02, confidence mínim de 0.4

rules <- apriori(transations, parameter = list(support = 0.02, confidence = 0.4))
## Apriori
## 
## Parameter specification:
##  confidence minval smax arem  aval originalSupport maxtime support minlen
##         0.4    0.1    1 none FALSE            TRUE       5    0.02      1
##  maxlen target   ext
##      10  rules FALSE
## 
## Algorithmic control:
##  filter tree heap memopt load sort verbose
##     0.1 TRUE TRUE  FALSE TRUE    2    TRUE
## 
## Absolute minimum support count: 336 
## 
## set item appearances ...[0 item(s)] done [0.00s].
## set transactions ...[3578 item(s), 16835 transaction(s)] done [0.21s].
## sorting and recoding items ... [198 item(s)] done [0.00s].
## creating transaction tree ... done [0.01s].
## checking subsets of size 1 2 3 done [0.00s].
## writing ... [51 rule(s)] done [0.00s].
## creating S4 object  ... done [0.00s].
inspect(head(sort(rules, by = "confidence"), 10))
##      lhs                                     rhs                                    support confidence     lift count
## [1]  {PINK REGENCY TEACUP AND SAUCER,                                                                                
##       ROSES REGENCY TEACUP AND SAUCER }   => {GREEN REGENCY TEACUP AND SAUCER}   0.02197802  0.8915663 22.40227   370
## [2]  {GREEN REGENCY TEACUP AND SAUCER,                                                                               
##       PINK REGENCY TEACUP AND SAUCER}     => {ROSES REGENCY TEACUP AND SAUCER }  0.02197802  0.8371041 18.84044   370
## [3]  {PINK REGENCY TEACUP AND SAUCER}     => {GREEN REGENCY TEACUP AND SAUCER}   0.02625483  0.8215613 20.64326   442
## [4]  {GREEN REGENCY TEACUP AND SAUCER}    => {ROSES REGENCY TEACUP AND SAUCER }  0.03070983  0.7716418 17.36710   517
## [5]  {PINK REGENCY TEACUP AND SAUCER}     => {ROSES REGENCY TEACUP AND SAUCER }  0.02465102  0.7713755 17.36110   415
## [6]  {GARDENERS KNEELING PAD CUP OF TEA } => {GARDENERS KNEELING PAD KEEP CALM } 0.02583903  0.7213930 17.12927   435
## [7]  {GREEN REGENCY TEACUP AND SAUCER,                                                                               
##       ROSES REGENCY TEACUP AND SAUCER }   => {PINK REGENCY TEACUP AND SAUCER}    0.02197802  0.7156673 22.39453   370
## [8]  {ROSES REGENCY TEACUP AND SAUCER }   => {GREEN REGENCY TEACUP AND SAUCER}   0.03070983  0.6911765 17.36710   517
## [9]  {DOLLY GIRL LUNCH BOX}               => {SPACEBOY LUNCH BOX }               0.02257202  0.6713781 17.82752   380
## [10] {ALARM CLOCK BAKELIKE GREEN}         => {ALARM CLOCK BAKELIKE RED }         0.03059103  0.6679637 13.21406   515

Amb aquesta ordenació per confidence, tenim una probabilitat alta que si compren {PINK REGENCY TEACUP AND SAUCER, ROSES REGENCY TEACUP AND SAUCER} compraran {GREEN REGENCY TEACUP AND SAUCER}. És un exemple, el mateix passa amb la resta de les 9 regles mostrades. En tots els casos el lift és molt alt, per tant no hi ha cap tipus d’aleatorietat en els resultats.

inspect(head(sort(rules, by = "support"), 10))
##      lhs                                   rhs                                   support confidence      lift count
## [1]  {GREEN REGENCY TEACUP AND SAUCER}  => {ROSES REGENCY TEACUP AND SAUCER } 0.03070983  0.7716418 17.367098   517
## [2]  {ROSES REGENCY TEACUP AND SAUCER } => {GREEN REGENCY TEACUP AND SAUCER}  0.03070983  0.6911765 17.367098   517
## [3]  {ALARM CLOCK BAKELIKE GREEN}       => {ALARM CLOCK BAKELIKE RED }        0.03059103  0.6679637 13.214064   515
## [4]  {ALARM CLOCK BAKELIKE RED }        => {ALARM CLOCK BAKELIKE GREEN}       0.03059103  0.6051704 13.214064   515
## [5]  {LUNCH BAG PINK POLKADOT}          => {LUNCH BAG RED RETROSPOT}          0.02863083  0.5446328  7.770248   482
## [6]  {LUNCH BAG RED RETROSPOT}          => {LUNCH BAG PINK POLKADOT}          0.02863083  0.4084746  7.770248   482
## [7]  {LUNCH BAG  BLACK SKULL.}          => {LUNCH BAG RED RETROSPOT}          0.02827443  0.4783920  6.825194   476
## [8]  {LUNCH BAG RED RETROSPOT}          => {LUNCH BAG  BLACK SKULL.}          0.02827443  0.4033898  6.825194   476
## [9]  {JUMBO BAG PINK POLKADOT}          => {JUMBO BAG RED RETROSPOT}          0.02684883  0.5955204  7.641453   452
## [10] {PINK REGENCY TEACUP AND SAUCER}   => {GREEN REGENCY TEACUP AND SAUCER}  0.02625483  0.8215613 20.643261   442

Si ordenem per support, veiem les regles que més cops s’han produït a tot el conjunt. La regla comprar {GREEN REGENCY TEACUP AND SAUCER} => {ROSES REGENCY TEACUP AND SAUCER } s’ha produït més de 517 cops. Aquestes regles tenen a més una probabilitat condicionada (confidence) relativament alta > 40%, algunes amb més d’un 80%.

El que queda força clar que els productes són similars i d’aquí el patró de compra.

Ara executarem el mateix algoritme baixant el support i pujan el confidence.

rules2 <- apriori(transations, parameter = list(support = 0.005, confidence = 0.9))
## Apriori
## 
## Parameter specification:
##  confidence minval smax arem  aval originalSupport maxtime support minlen
##         0.9    0.1    1 none FALSE            TRUE       5   0.005      1
##  maxlen target   ext
##      10  rules FALSE
## 
## Algorithmic control:
##  filter tree heap memopt load sort verbose
##     0.1 TRUE TRUE  FALSE TRUE    2    TRUE
## 
## Absolute minimum support count: 84 
## 
## set item appearances ...[0 item(s)] done [0.00s].
## set transactions ...[3578 item(s), 16835 transaction(s)] done [0.30s].
## sorting and recoding items ... [1161 item(s)] done [0.01s].
## creating transaction tree ... done [0.01s].
## checking subsets of size 1 2 3 4 5 6 done [0.07s].
## writing ... [221 rule(s)] done [0.00s].
## creating S4 object  ... done [0.01s].
inspect(head(sort(rules2, by = "confidence"), 10))
##      lhs                                   rhs                            support confidence     lift count
## [1]  {HERB MARKER BASIL,                                                                                   
##       HERB MARKER CHIVES ,                                                                                 
##       HERB MARKER MINT,                                                                                    
##       HERB MARKER ROSEMARY,                                                                                
##       HERB MARKER THYME}                => {HERB MARKER PARSLEY}      0.007722008  0.9923664 90.79613   130
## [2]  {GREEN REGENCY TEACUP AND SAUCER,                                                                     
##       REGENCY TEA PLATE PINK,                                                                              
##       REGENCY TEA PLATE ROSES }         => {REGENCY TEA PLATE GREEN } 0.005524206  0.9893617 70.57587    93
## [3]  {GREEN REGENCY TEACUP AND SAUCER,                                                                     
##       REGENCY TEA PLATE PINK,                                                                              
##       REGENCY TEA PLATE ROSES ,                                                                            
##       ROSES REGENCY TEACUP AND SAUCER } => {REGENCY TEA PLATE GREEN } 0.005227205  0.9887640 70.53323    88
## [4]  {HERB MARKER CHIVES ,                                                                                 
##       HERB MARKER MINT,                                                                                    
##       HERB MARKER ROSEMARY,                                                                                
##       HERB MARKER THYME}                => {HERB MARKER PARSLEY}      0.008197208  0.9857143 90.18750   138
## [5]  {HERB MARKER BASIL,                                                                                   
##       HERB MARKER CHIVES ,                                                                                 
##       HERB MARKER ROSEMARY,                                                                                
##       HERB MARKER THYME}                => {HERB MARKER PARSLEY}      0.007959608  0.9852941 90.14906   134
## [6]  {HERB MARKER BASIL,                                                                                   
##       HERB MARKER CHIVES ,                                                                                 
##       HERB MARKER PARSLEY,                                                                                 
##       HERB MARKER ROSEMARY}             => {HERB MARKER THYME}        0.007959608  0.9852941 91.13971   134
## [7]  {HERB MARKER BASIL,                                                                                   
##       HERB MARKER CHIVES ,                                                                                 
##       HERB MARKER MINT,                                                                                    
##       HERB MARKER THYME}                => {HERB MARKER PARSLEY}      0.007840808  0.9850746 90.12897   132
## [8]  {HERB MARKER BASIL,                                                                                   
##       HERB MARKER CHIVES ,                                                                                 
##       HERB MARKER MINT,                                                                                    
##       HERB MARKER ROSEMARY}             => {HERB MARKER PARSLEY}      0.007840808  0.9850746 90.12897   132
## [9]  {HERB MARKER BASIL,                                                                                   
##       HERB MARKER CHIVES ,                                                                                 
##       HERB MARKER MINT,                                                                                    
##       HERB MARKER PARSLEY,                                                                                 
##       HERB MARKER THYME}                => {HERB MARKER ROSEMARY}     0.007722008  0.9848485 88.66270   130
## [10] {HERB MARKER BASIL,                                                                                   
##       HERB MARKER CHIVES ,                                                                                 
##       HERB MARKER MINT,                                                                                    
##       HERB MARKER PARSLEY,                                                                                 
##       HERB MARKER ROSEMARY}             => {HERB MARKER THYME}        0.007722008  0.9848485 91.09848   130

Ara hem modificat l’algoritme per a produir regles amb un confidence molt alt, encara que no es produeixin gaires cops en el total del dataset(support més baix). Tenim similaritat de productes, però encara així seria una bona forma de treure campanyes de marketing.

inspect(head(sort(rules2, by = "support"), 10))
##      lhs                                rhs                            support confidence     lift count
## [1]  {POPPY'S PLAYHOUSE BEDROOM ,                                                                       
##       POPPY'S PLAYHOUSE LIVINGROOM } => {POPPY'S PLAYHOUSE KITCHEN} 0.01069201  0.9000000 44.56324   180
## [2]  {HERB MARKER THYME}             => {HERB MARKER ROSEMARY}      0.01021681  0.9450549 85.08021   172
## [3]  {HERB MARKER ROSEMARY}          => {HERB MARKER THYME}         0.01021681  0.9197861 85.08021   172
## [4]  {HERB MARKER PARSLEY}           => {HERB MARKER ROSEMARY}      0.01015741  0.9293478 83.66615   171
## [5]  {HERB MARKER ROSEMARY}          => {HERB MARKER PARSLEY}       0.01015741  0.9144385 83.66615   171
## [6]  {HERB MARKER MINT}              => {HERB MARKER ROSEMARY}      0.01009801  0.9042553 81.40716   170
## [7]  {HERB MARKER ROSEMARY}          => {HERB MARKER MINT}          0.01009801  0.9090909 81.40716   170
## [8]  {HERB MARKER BASIL}             => {HERB MARKER ROSEMARY}      0.01003861  0.9184783 82.68760   169
## [9]  {HERB MARKER ROSEMARY}          => {HERB MARKER BASIL}         0.01003861  0.9037433 82.68760   169
## [10] {HERB MARKER PARSLEY}           => {HERB MARKER MINT}          0.00997921  0.9130435 81.76110   168

El mateix ordenat per support. Els lifts són altíssims, per sobre d’1 ens donaríem per satisfets.

inspect(head(sort(rules2, by = "lift"), 10))
##      lhs                            rhs                             support confidence      lift count
## [1]  {DOLLY GIRL CHILDRENS CUP,                                                                       
##       SPACEBOY CHILDRENS BOWL}   => {DOLLY GIRL CHILDRENS BOWL} 0.005405405  0.9680851 120.72380    91
## [2]  {DOLLY GIRL CHILDRENS BOWL,                                                                      
##       SPACEBOY CHILDRENS CUP}    => {DOLLY GIRL CHILDRENS CUP}  0.005049005  0.9550562 116.50993    85
## [3]  {DOLLY GIRL CHILDRENS BOWL,                                                                      
##       SPACEBOY CHILDRENS CUP}    => {SPACEBOY CHILDRENS BOWL}   0.005049005  0.9550562  98.64031    85
## [4]  {HERB MARKER BASIL,                                                                              
##       HERB MARKER MINT,                                                                               
##       HERB MARKER PARSLEY,                                                                            
##       HERB MARKER ROSEMARY,                                                                           
##       HERB MARKER THYME}         => {HERB MARKER CHIVES }       0.007722008  0.9219858  95.81254   130
## [5]  {HERB MARKER MINT,                                                                               
##       HERB MARKER PARSLEY,                                                                            
##       HERB MARKER ROSEMARY,                                                                           
##       HERB MARKER THYME}         => {HERB MARKER CHIVES }       0.008197208  0.9139073  94.97302   138
## [6]  {HERB MARKER BASIL,                                                                              
##       HERB MARKER MINT,                                                                               
##       HERB MARKER PARSLEY,                                                                            
##       HERB MARKER THYME}         => {HERB MARKER CHIVES }       0.007840808  0.9103448  94.60281   132
## [7]  {HERB MARKER BASIL,                                                                              
##       HERB MARKER PARSLEY,                                                                            
##       HERB MARKER ROSEMARY,                                                                           
##       HERB MARKER THYME}         => {HERB MARKER CHIVES }       0.007959608  0.9054054  94.08951   134
## [8]  {HERB MARKER BASIL,                                                                              
##       HERB MARKER MINT,                                                                               
##       HERB MARKER PARSLEY,                                                                            
##       HERB MARKER ROSEMARY}      => {HERB MARKER CHIVES }       0.007840808  0.9041096  93.95485   132
## [9]  {HERB MARKER MINT,                                                                               
##       HERB MARKER PARSLEY,                                                                            
##       HERB MARKER THYME}         => {HERB MARKER CHIVES }       0.008375408  0.9038462  93.92747   141
## [10] {HERB MARKER BASIL,                                                                              
##       HERB MARKER CHIVES ,                                                                            
##       HERB MARKER PARSLEY,                                                                            
##       HERB MARKER ROSEMARY}      => {HERB MARKER THYME}         0.007959608  0.9852941  91.13971   134

La conclusió que podem treure d’aquest estudi és que a partir de l’algoritme d’associació podem veure patrons de conducta en les compres que es repeteixen i podem trobar associacions amb una fiabilitat molt alta.


3 Clustering (concepte de distància): K-means o centroides


Aquest algoritme primerament fixa un nombre de clusters (o centres de clusters) i a partir d’aquest nombre construeix els grups. Per a executar el mètode de k-means o centroides partim de la base que no coneixem el nombre òptim de clústers. Provarem amb diversos centres (del 2 al 10) per veure quina aproximació és millor.

Com funciona Kmeans?

Havíem preparat el model per fer dos estudis: de productes com de clients. Començarem pels clients.

3.1 Estudi de customers

customer.data <- read.csv('retail-customer.csv', sep = ',')
summary(customer.data)
##   NumInvoices     NumCancellations    TotalCost          NumStocks      
##  Min.   :  0.00   Min.   : 0.0000   Min.   : -4287.6   Min.   :   0.00  
##  1st Qu.:  1.00   1st Qu.: 0.0000   1st Qu.:   293.4   1st Qu.:  15.00  
##  Median :  2.00   Median : 0.0000   Median :   648.1   Median :  35.00  
##  Mean   :  4.24   Mean   : 0.8358   Mean   :  1898.5   Mean   :  61.03  
##  3rd Qu.:  5.00   3rd Qu.: 1.0000   3rd Qu.:  1611.7   3rd Qu.:  77.00  
##  Max.   :210.00   Max.   :47.0000   Max.   :279489.0   Max.   :1787.00  
##       Days         DiffInvoices    
##  Min.   :  0.00   Min.   : -6.000  
##  1st Qu.: 16.10   1st Qu.:  1.000  
##  Median : 49.90   Median :  2.000  
##  Mean   : 91.59   Mean   :  3.404  
##  3rd Qu.:142.93   3rd Qu.:  4.000  
##  Max.   :373.10   Max.   :196.000

Anem a fer l’estudi amb 4 variables:

  • Days: Dies des de la darrera compra (respecte a la data del darrer producte comprat)
  • NumInvoices: Nombre de factures o compres realitzades
  • TotalCost: Cost gastat total
  • NumStocks: Nombre de productes comprats en total
# filtrem els camps a treballar
customer.data.clean <- customer.data[,c('Days','NumInvoices','TotalCost','NumStocks')]
summary(customer.data.clean)
##       Days         NumInvoices       TotalCost          NumStocks      
##  Min.   :  0.00   Min.   :  0.00   Min.   : -4287.6   Min.   :   0.00  
##  1st Qu.: 16.10   1st Qu.:  1.00   1st Qu.:   293.4   1st Qu.:  15.00  
##  Median : 49.90   Median :  2.00   Median :   648.1   Median :  35.00  
##  Mean   : 91.59   Mean   :  4.24   Mean   :  1898.5   Mean   :  61.03  
##  3rd Qu.:142.93   3rd Qu.:  5.00   3rd Qu.:  1611.7   3rd Qu.:  77.00  
##  Max.   :373.10   Max.   :210.00   Max.   :279489.0   Max.   :1787.00

Eliminarem els outliers de totalCost, ja que crec que poden desvirtuar la realitat. Per a trobar valors extrems anem a aplicar la idea dels IQR (interquartile ranges):

Per a una determinada variable contínua, els outliers són aquelles observacions que es troben fora de 1.5 * IQR, on IQR, el “Inter Quartile Range” és la diferència entre el Q3 i el Q1:

  • Interquartile range, IQR = Q3 - Q1
  • lower = Q1 - 1.5 * IQR
  • Upper = Q3 + 1.5 * IQR
# funció per treure els límits dels valors atípics:
outliersLimits <- function(x) {
    limits <- c("above", "under")
    limits$above <- quantile(x, 0.75, type=6) + 1.5 * (quantile(x, 0.75, type=6) - quantile(x, 0.25, type=6))
    limits$under <- quantile(x, 0.25, type=6) - 1.5 * (quantile(x, 0.75, type=6) - quantile(x, 0.25, type=6))
  return(limits)
}

limits <- outliersLimits(customer.data.clean$TotalCost)

# treiem aquest valors atípics de total cost. Per el límit superiot anem a agafar aquells que tenen una despesa superior a 0
customer.data.clean <-subset(customer.data.clean, TotalCost > 0)
customer.data.clean <-subset(customer.data.clean, TotalCost <= limits$above)
summary(customer.data.clean)
##       Days         NumInvoices       TotalCost        NumStocks     
##  Min.   :  0.00   Min.   : 1.000   Min.   :   2.9   Min.   :  1.00  
##  1st Qu.: 20.80   1st Qu.: 1.000   1st Qu.: 275.3   1st Qu.: 15.00  
##  Median : 56.10   Median : 2.000   Median : 581.8   Median : 31.00  
##  Mean   : 97.27   Mean   : 2.923   Mean   : 869.3   Mean   : 48.75  
##  3rd Qu.:155.00   3rd Qu.: 4.000   3rd Qu.:1223.1   3rd Qu.: 65.00  
##  Max.   :373.10   Max.   :39.000   Max.   :3580.1   Max.   :580.00

Com la magnitud dels valors difereix notablement entre les variables, les escalem abans d’aplicar el clustering.

df1 <- scale(customer.data.clean)

Mostrem en una gràfica els valors de la silueta mitjana de cada prova per comprovar quin nombre de clusters és el millor. Com més alt és el valor de la silueta, millor

avg_sil <- function(k) {
  km.res <- kmeans(df1, centers = k, nstart = 25, iter.max=25)
  ss <- silhouette(km.res$cluster, dist(df1))
  mean(ss[, 3])
}

k.values <- 2:12

avg_sil_values <- map_dbl(k.values, avg_sil)

plot(k.values, avg_sil_values,
       type = "b", pch = 19, frame = FALSE, 
       xlab = "Nombre de clusters K",
       ylab = "Avg. Silhouettes")

avg_sil_values
##  [1] 0.4549749 0.4022606 0.3755009 0.3620711 0.3456050 0.3495608 0.3122243
##  [8] 0.2937851 0.2958048 0.2938006 0.2971479

El millor valor obtingut és de 2 clusters o tipus de clients, però del 3 al 7 estan relativament molt a prop quant a valor de Avg. Silhouettes.

Podem executar el mateix procés amb la funció fviz_nbclust

fviz_nbclust(df1, kmeans, method = "silhouette")

Res de nou. Surten 2 clusters. I del 3 al 7, bons valors també.

Un altra forma d’avaluar quin és el millor nombre de clústers és considerar el millor model amb aquestes condicions:

  • Ofereix la menor suma dels quadrats de les distàncies dels punts de cada grup respecte al seu centre (withinss). Això vol dir que els elements dels grups estan junts.

  • La separació més gran entre centres de grups (betweenss). Els grups están separats.

És una idea conceptualment similar a la silueta. Una manera comú de fer la selecció del nombre de clústers consisteix a aplicar el mètode elbow (colze). Es fa una selecció del nombre de clústers a partir de la inspecció de la gràfica que s’obté l’iterar amb el mateix conjunt de dades per a distints valors del nombre de clústers. Es seleccionarà el valor que es troba en el colze de la corba.

wss <- function(k) {
  kmeans(df1, k, nstart = 25, iter.max = 25 )$tot.withinss
}

# fem fins a 15 iteracions
k.values <- 1:15


wss_values <- map_dbl(k.values, wss)

plot(k.values, wss_values,
       type="b", pch = 19, frame = FALSE, 
       xlab="Nombre de clusters K",
       ylab="Total within-clusters sum of squares")

# la funció fviz_nbclust
fviz_nbclust(df1, kmeans, method = "wss")

A partir de la corba obtinguda podem veure com a mesura que s’augmenta la quantitat de centroides, el valor de “tot.tot.withinss” disminueix, ja que la separació entre elements dels clusters és inferior. La idea és trobar un “colze”. El colze es troba on ja no es produeixen variacions importants en augmentar ‘Nombre de clusters’. El valor que jo dedueixo és el 3, és on tenim el colze.

També es pot fer servir la funció kmeansruns del paquet fpc que executarà l’algoritme kmeans com un conjunt de valors i selecciona el valor del nombre de clústers que millor funcioni d’acord amb dos criteris: la silueta mitja (asw) i Calinski-Harabasz (“ch”).

fit_ch  <- kmeansruns(df1, krange = 1:10, criterion = "ch", iter.max=25, runs=25)
plot(1:10,fit_ch$crit,type="o",col="blue",pch=0,xlab="Nombre de clústers",ylab="Criteri Calinski-Harabasz")

Amb aquest criteri estem obtenint 3 clusters també.

fit_asw <- kmeansruns(df1, krange = 1:10, criterion = "asw", iter.max=25, runs=25)
plot(1:10,fit_asw$crit,type="o",col="blue",pch=0,xlab="Nombre de clústers",ylab="Criteri silueta mitja")

Realitzem comparacions visuals amb el nombre de clusters que hem trobat en la majoria d’estudis (2 i 3 clusters)

# realitzem una divisió en 2 clusters
clusters <- kmeans(df1, centers = 2, iter.max = 25, nstart = 25)
clusters$centers
##         Days NumInvoices  TotalCost  NumStocks
## 1 -0.6280280   1.2333845  1.3431643  1.1268210
## 2  0.2151878  -0.4226074 -0.4602225 -0.3860945
fviz_cluster(clusters, data = df1, geom = "point", repel = TRUE)

El mapa ens diferencia 2 grups de clients força clars.

fviz_cluster(clusters, data = df1, geom = "point", choose.vars = c("Days", "NumInvoices"), repel = TRUE, ggtheme = theme_minimal())

Si mostrem el nombre de factures x dies. Veiem una sèrie de clients potencialment molt bons, on han comprat relativament fa poc i tenen moltes factures. Són clients habituals.

fviz_cluster(clusters, data = df1, geom = "point", choose.vars = c("Days", "TotalCost"), repel = TRUE, ggtheme = theme_minimal())

En aquesta mostra tenim els millors clients en el grupuscle superior-esquerra. Més despesa i menys dies des de la darrera compra.

fviz_cluster(clusters, data = df1, geom = "point", choose.vars = c("NumInvoices", "TotalCost"), repel = TRUE, ggtheme = theme_minimal())

I en l’estudi de les dues variables, nombre de factures x cost. Diferenciem clients que gasten molt amb poques factures (part superior - esquerra) i altres que gasten força amb més factures (part superior - dreta).

Fem l’estudi visual amb 3 clusters:

# realitzem una divisió en 3 clusters
clusters <- kmeans(df1, centers = 3, iter.max = 25, nstart = 25)
clusters$centers
##         Days NumInvoices  TotalCost  NumStocks
## 1 -0.4628347  -0.2483648 -0.2928991 -0.2352042
## 2  1.5221179  -0.5952659 -0.6190990 -0.5331632
## 3 -0.6343581   1.4406513  1.5944061  1.3257665
fviz_cluster(clusters, data = df1, geom = "point", repel = TRUE)

Un dels dos clusters que teníem s’ha dividit en dos (cluster 1 i 2).

set.seed(123)
fviz_cluster(clusters, data = df1, geom = "point", choose.vars = c("Days", "NumInvoices"), repel = TRUE, ggtheme = theme_minimal())

Mostrem el nombre de factures x dies:

  • el grup 1, fa relativament poc que han realitzat compra però amb menys factures.
  • el grup 2, poques compres i fa més dies que no compren.
  • el grup 3, és potencialment el millor, tenen diverses factures i fa pocs dies que han comprat, reincideixen.
set.seed(123)
fviz_cluster(clusters, data = df1, geom = "point", choose.vars = c("Days", "TotalCost"), repel = TRUE, ggtheme = theme_minimal())

En aquesta mostra tenim, Days vs Cost:

  • el grup 1, fa poc que han realitzat compra però amb menys cost
  • el grup 2, fa més temps que no compren i han gastat menys.
  • els millors clients en el grup 3.
set.seed(123)
fviz_cluster(clusters, data = df1, geom = "point", choose.vars = c("NumInvoices", "TotalCost"), repel = TRUE, ggtheme = theme_minimal())

I en l’estudi de les dues variables, nombre de factures x cost:

  • hi ha dos grups que pràcticament en superposen.
  • un tercer que diferencia el nombre de factures (o compres) que han realitzat.

En l’estudi hem optat per a realitzar les gràfiques visuals amb 2 i 3 clusters. De tota manera, els primers algoritmes ens estaven donant possibilitat de més agrupacions fins a un total de 7.

3.2 Estudi de productes

Ara farem un estudi similar amb els productes.

retail.product.data <- read.csv('retail-product.csv', sep = ',')
summary(retail.product.data)
##    UnitPrice       TotalStockUnits        Days        TotalStockPrice    
##  Min.   :  0.000   Min.   :-1460.0   Min.   :  0.00   Min.   : -86541.1  
##  1st Qu.:  1.010   1st Qu.:   67.0   1st Qu.:  0.90   1st Qu.:    142.7  
##  Median :  1.950   Median :  414.5   Median :  3.10   Median :    750.7  
##  Mean   :  3.818   Mean   : 1395.3   Mean   : 46.06   Mean   :   2861.7  
##  3rd Qu.:  3.880   3rd Qu.: 1437.0   3rd Qu.: 32.02   3rd Qu.:   2373.4  
##  Max.   :744.150   Max.   :53215.0   Max.   :373.10   Max.   :1064825.1

Anem a fer l’estudi amb 3 variables:

  • UnitPrice: Preu unitari del producte.
  • Days: Dies des de la darrera compra (respecte a la data del darrer producte comprat)
  • TotalStockUnits: Unitats comprades totals
product.data.clean <- retail.product.data[,c('UnitPrice', 'Days','TotalStockUnits')]
summary(product.data.clean)
##    UnitPrice            Days        TotalStockUnits  
##  Min.   :  0.000   Min.   :  0.00   Min.   :-1460.0  
##  1st Qu.:  1.010   1st Qu.:  0.90   1st Qu.:   67.0  
##  Median :  1.950   Median :  3.10   Median :  414.5  
##  Mean   :  3.818   Mean   : 46.06   Mean   : 1395.3  
##  3rd Qu.:  3.880   3rd Qu.: 32.02   3rd Qu.: 1437.0  
##  Max.   :744.150   Max.   :373.10   Max.   :53215.0
# eliminem valors negatius o igual a 0
product.data.clean <- subset(product.data.clean, product.data.clean$TotalStockUnits > 0)
product.data.clean <- subset(product.data.clean, product.data.clean$UnitPrice > 0)
product.data.clean <- subset(product.data.clean, product.data.clean$Days > 0)
summary(product.data.clean)
##    UnitPrice            Days        TotalStockUnits
##  Min.   :  0.040   Min.   :  0.10   Min.   :    1  
##  1st Qu.:  0.990   1st Qu.:  1.00   1st Qu.:   65  
##  Median :  1.950   Median :  3.80   Median :  371  
##  Mean   :  3.699   Mean   : 46.63   Mean   : 1227  
##  3rd Qu.:  3.840   3rd Qu.: 34.90   3rd Qu.: 1264  
##  Max.   :744.150   Max.   :373.10   Max.   :53215

Eliminarem els outliers de total cost, ja que crec que poden desvirtuar la realitat. Funció per treure els límits dels valors atípics:

outliersLimits <- function(x) {
    limits <- c("above", "under")
    limits$above <- quantile(x, 0.75, type=6) + 1.5 * (quantile(x, 0.75, type=6) - quantile(x, 0.25, type=6))
    limits$under <- quantile(x, 0.25, type=6) - 1.5 * (quantile(x, 0.75, type=6) - quantile(x, 0.25, type=6))
  return(limits)
}

limits <- outliersLimits(product.data.clean$UnitPrice)

# treiem aquest valors atípics de total cost
product.data.clean <-subset(product.data.clean, UnitPrice <= limits$above)

summary(product.data.clean)
##    UnitPrice          Days        TotalStockUnits
##  Min.   :0.040   Min.   :  0.10   Min.   :    1  
##  1st Qu.:0.870   1st Qu.:  1.00   1st Qu.:   75  
##  Median :1.660   Median :  3.80   Median :  426  
##  Mean   :2.327   Mean   : 46.23   Mean   : 1306  
##  3rd Qu.:3.150   3rd Qu.: 32.90   3rd Qu.: 1375  
##  Max.   :8.110   Max.   :373.10   Max.   :53215

Com que no volem que l’algoritme de clúster no depengui d’una unitat variable arbitrària, comencem escalant / estandarditzant les dades mitjançant la funció scale.

df2 <- scale(product.data.clean)

Anem a fer servir el mètode de la silueta per a trobar el millor nombre de clusters. En resum, l’enfocament mitjà de la silueta mesura la qualitat d’una agrupació. És a dir, determina el bé que cada objecte es troba dins del seu clúster. Una gran amplada mitjana de silueta indica una bona agrupació. El mètode de silueta mitjana calcula la silueta mitjana d’observacions per a diferents valors de k.

avg_sil <- function(k) {
  km.res <- kmeans(df2, centers = k, nstart = 25, iter.max=25)
  ss <- silhouette(km.res$cluster, dist(df2))
  mean(ss[, 3])
}

k.values <- 2:15

avg_sil_values <- map_dbl(k.values, avg_sil)

plot(k.values, avg_sil_values,
       type = "b", pch = 19, frame = FALSE, 
       xlab = "Nombre de clusters K",
       ylab = "Avg. Silhouettes")

Segons el mètode de la silueta obtenim el millor nombre de clusters entre 4 i 5 tipus de productes diferents.

Podem executar el mateix procés amb la funció fviz_nbclust

fviz_nbclust(df2, kmeans, method = "silhouette")

4 és el nombre de clusters, seguit del 5 altre cop.

# function to compute total within-cluster sum of square 
wss <- function(k) {
  kmeans(df2, k, nstart = 25,iter.max=25)$tot.withinss
}

# Compute and plot wss for k = 1 to k = 15
k.values <- 1:15

# extract wss for 2-15 clusters
wss_values <- map_dbl(k.values, wss)

plot(k.values, wss_values,
       type="b", pch = 19, frame = FALSE, 
       xlab="Number of clusters K",
       ylab="Total within-clusters sum of squares")

Tenim el colze on ja no es produeixen variacions importants en augmentar ‘Nombre de clusters’. El valor que jo dedueixo és entre 4 i 5, altre cop.

Afortunadament, aquest procés per calcular el “mètode del colze” i la silueta s’ha agrupat en una única funció (fviz_nbclust):

fviz_nbclust(df2, kmeans, method = "wss")

No hi ha molt dubte. Tenim entre 4 tipus de productes diferents segons dies des de la darrera compra, unitats comprades i cost del producte.

Anem a fer servir la funció kmeansruns del paquet fpc que executarà l’algoritme kmeans com un conjunt de valors i selecciona el valor del nombre de clústers que millor funcioni d’acord amb dos criteris: la silueta mitja (asw) i Calinski-Harabasz (“ch”).

fit_asw <- kmeansruns(df2, krange = 2:20, criterion = "asw", iter.max=25, runs=10)
plot(1:20,fit_asw$crit,type="o",col="blue",pch=0,xlab="Nombre de clústers",ylab="Criteri silueta mitja")

Tornem a tenir 4 o 5 clusters. Farem estudis visuals amb aquests valors, ja que els resultats són molt més clars que amb els customers.

# realitzem una divisió en 4 clusters
clusters4 <- kmeans(df2, centers = 4, iter.max=25, nstart = 25)
clusters4$centers
##     UnitPrice       Days TotalStockUnits
## 1 -0.56788319 -0.5063242      3.62724749
## 2  1.51580309 -0.2952085     -0.27705397
## 3 -0.47872208 -0.3864604     -0.05177185
## 4 -0.03026779  2.2455655     -0.42090266
fviz_cluster(clusters4, data = df2, geom = "point", repel = TRUE)

set.seed(123)
fviz_cluster(clusters4, data = df2, geom = "point", choose.vars = c("UnitPrice", "TotalStockUnits"), repel = TRUE, ggtheme = theme_minimal())

Unitats venudes totals x Preu unitari:

  • Grup 1: es venen força unitats i tenim preus baixos i mitjans.
  • Grup 2: tenim preus mitjans i alts, però es venen poques unitats.
  • Grup 3: preus mitjans i baixos i poques unitats venudes.
  • Grup 4 és un grup molt solapat amb els altres. Segurament la diferencia es troba en l’altre variable no inclosa (Days).
set.seed(123)
fviz_cluster(clusters4, data = df2, geom = "point", choose.vars = c("UnitPrice", "Days"), repel = TRUE, ggtheme = theme_minimal())

Preu unitari x dies des de la darrera compra:

  • Grup 1: es troba molt encavalcat amb el 3 i 4.
  • Grup 2: fa poc que s’han comprat i tenen preus mitjans-alts.
  • Grup 3: fa poc que s’han comprat i tenen preus mitjans-baixos.
  • Grup 4: fa dies que no es compra i preus en tot el ventall.
set.seed(123)
fviz_cluster(clusters4, data = df2, geom = "point", choose.vars = c("TotalStockUnits", "Days"), repel = TRUE, ggtheme = theme_minimal())

Dies des de la darrera compra x Unitats venudes:

  • Grup 1: s’han venut relativament de mig a moltes unitats i fa poc de la darrera compra.
  • Grup 2: s’han venut poques unitats i la darrera compra és recent.
  • Grup 3: superposat amb els altres. Depèn de l’altre variable.
  • Grup 4: s’han venut poques unitats i fa temps.
# realitzem una divisió en 5 clusters
clusters5 <- kmeans(df2, centers = 5, iter.max=25, nstart = 25)
clusters5$centers
##     UnitPrice       Days TotalStockUnits
## 1 -0.50342869 -0.5161059       9.7281381
## 2 -0.47372940 -0.3825041      -0.1214569
## 3  1.51933002 -0.2939350      -0.2847242
## 4 -0.55164654 -0.4897042       2.2492208
## 5 -0.03026779  2.2455655      -0.4209027
summary(df2)
##    UnitPrice            Days         TotalStockUnits   
##  Min.   :-1.2337   Min.   :-0.5265   Min.   :-0.46919  
##  1st Qu.:-0.7859   1st Qu.:-0.5163   1st Qu.:-0.44259  
##  Median :-0.3597   Median :-0.4843   Median :-0.31645  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.00000  
##  3rd Qu.: 0.4442   3rd Qu.:-0.1521   3rd Qu.: 0.02471  
##  Max.   : 3.1201   Max.   : 3.7312   Max.   :18.65567
fviz_cluster(clusters5, data = df2, geom = "point", repel = TRUE)

Realment no ens aporta gaire generar un altre cluster. Amb 4 clusters ens quedem.


4 Clustering (concepte de distància): Mètode agregació jerarquica incremental (Hierarchical Agglomerative Clustering - HAC)


Funciona amb l’algoritme dels veïns més propers (k-nearest neighbours) de la següent forma:

De formes de mesurar les distàncies entre grups tenim de diversos tipus. Provarem només: “average”, “complete” i “ward”.

4.1 Estudi de customers

Fem servir el dataframe: ‘customer.data.clean’ que hem generat a Clustering amb K-means o centroides:

summary(customer.data.clean)
##       Days         NumInvoices       TotalCost        NumStocks     
##  Min.   :  0.00   Min.   : 1.000   Min.   :   2.9   Min.   :  1.00  
##  1st Qu.: 20.80   1st Qu.: 1.000   1st Qu.: 275.3   1st Qu.: 15.00  
##  Median : 56.10   Median : 2.000   Median : 581.8   Median : 31.00  
##  Mean   : 97.27   Mean   : 2.923   Mean   : 869.3   Mean   : 48.75  
##  3rd Qu.:155.00   3rd Qu.: 4.000   3rd Qu.:1223.1   3rd Qu.: 65.00  
##  Max.   :373.10   Max.   :39.000   Max.   :3580.1   Max.   :580.00

Abans de tot, normalitzem les dades i calculem la distància euclidiana entre els elements.

## creem la funció per normalitzar
ni <-function (x) {(x -min (x)) / (max (x) -min (x))} 
customer.data.norm <- as.data.frame(lapply(customer.data.clean[,1:4], ni))
# primer fem la matriu de la distància euclideana dels elements de la taula de llavors.
customer.data.euc <- dist(customer.data.norm, method = "euclidean")

4.1.1 Enllaç complet

Estratègia de la distància màxima o similitud mínima. En aquest mètode, també conegut (complete linkage), es considera que la distància o similitud entre dos clústers cal mesurar-la atenent a elements més dispars, és a dir, la distància o similitud entre clústers ve dada, respectivament, per la màxima distància (o mínima similitud) entre components dels clústers.

# fem servir hclust per a realitzar les agrupacions
hc1 <- hclust(customer.data.euc, method = "complete")
plot(hc1, 
     main = "Complete linkage", 
     xlab = "Customers",
     ylab = "",
     cex = 0.04,
     sub = "")
abline(h = 0.78, col = "red")
abline(h = 1.1, col = "blue")

El dendograma ens està mostrant les distàncies euclidianes entre els diferents elements (customers) de la mostra i la separació gràfica conformes en van agrupant els clusters. Les interpretació és subjectiva però agafariem 4 i 7 clusters.

4.1.2 Enllaç promig

En aquesta estratègia la distància, o similitud, de l’clúster Ci amb el Cj s’obté com la mitjana aritmètica entre la distància, o similitud, dels elements d’aquests clústers.

hc2 <- hclust(customer.data.euc, method = "average")

plot(hc2, 
     main = "Average linkage", 
     xlab = "Customers",
     ylab = "",
     cex = 0.1,
     sub = "")
abline(h = 0.4, col = "red")
abline(h = 0.55, col = "blue")

Similar, entre 3 i 7 grups de customers.

4.1.3 Mètode Ward

El mètode de Ward apunta a minimitzar la variància total dins de el grup. A cada pas, es fusionen el parell de clústers amb una distància mínima entre els clústers. En altres paraules, forma grups d’una manera que minimitza la pèrdua associada amb cada grup. Té molt bona efectivitat a l’hora d’agrupar.

hc3 <- hclust(customer.data.euc, method = "ward.D")

plot(hc3, 
     main = "Ward Method", 
     xlab = "Customers",
     ylab = "",
     cex = 0.1,
     sub = "")
abline(h = 100, col = "red")
abline(h = 55, col = "blue")

Amb el mètode ward veiem més clarament 4 o 5 grups de clients. Els altres clusters es troben molt més allunyats en l’agrupació. Ens quedariem amb 4 o 5 grups de customers.

Anem a reduir la quantitat d’elements de la mostra i apliquem algoritme ward

set.seed(123)
sample <- customer.data.norm %>% sample_frac(0.10)
sample.euc <- dist(sample, method = "euclidean")
hc3 <- hclust(sample.euc, method = "ward.D")

plot(hc3, 
     main = "Ward Method", 
     xlab = "Products",
     ylab = "",
     cex = 0.1,
     sub = "")
abline(h = 12, col = "red")
abline(h = 7, col = "green")

Amb menys mostra estem treient entre 3 i 5 clusters.

Dels 3 mètodes que hem fet servir el que dona una més bona interpretació és ward on podem veure entre 3 i 5 clusters.

4.2 Estudi de productes

Fem servir el dataframe: ‘product.data.clean’ que hem generat a Clustering amb K-means o centroides:

summary(product.data.clean)
##    UnitPrice          Days        TotalStockUnits
##  Min.   :0.040   Min.   :  0.10   Min.   :    1  
##  1st Qu.:0.870   1st Qu.:  1.00   1st Qu.:   75  
##  Median :1.660   Median :  3.80   Median :  426  
##  Mean   :2.327   Mean   : 46.23   Mean   : 1306  
##  3rd Qu.:3.150   3rd Qu.: 32.90   3rd Qu.: 1375  
##  Max.   :8.110   Max.   :373.10   Max.   :53215

Abans de tot, normalitzem les dades i calculem la distància euclidiana entre els elements.

## creem la funció per normalitzar
ni <-function (x) {(x -min (x)) / (max (x) -min (x))} 
product.data.norm <- as.data.frame(lapply(product.data.clean[,1:3], ni))
# primer fem la matriu de la distància euclideana dels elements de la taula de llavors.
product.data.euc <- dist(product.data.norm, method = "euclidean")

4.2.1 Mètode Ward

Executarem directament amb el mètode ward on veiem que tenim 4 clusters diferenciats.

hc3 <- hclust(product.data.euc, method = "ward.D")

plot(hc3, 
     main = "Ward Method", 
     xlab = "Products",
     ylab = "",
     cex = 0.1,
     sub = "")
abline(h = 70, col = "red")
abline(h = 50, col = "blue")

El dendograma ens està mostrant les distàncies euclidianes entre els diferents elements (productes) de la mostra. Mostra la separació dels elements i dels clusters que es van formant. Sembla distingir entre 4 i 5 clusters. La distància fins a 3 clusters ja és molt llunyana.

Reduïm la dimensionalitat de la mostra i apliquem algoritme ward.

set.seed(123)
sample <- product.data.norm %>% sample_frac(0.10)
sample.euc <- dist(sample, method = "euclidean")
hc3 <- hclust(sample.euc, method = "ward.D")

plot(hc3, 
     main = "Ward Method", 
     xlab = "Products",
     ylab = "",
     cex = 0.1,
     sub = "")
abline(h = 8, col = "red")
abline(h = 6, col = "green")

Podríem fer entre 4 i 5 clusters. Resultat molt similar al K-means o centroides.


5 Clustering (concepte probabilístic): Mètodes d’agregació probabilítics


Els enfocaments dels mètodes probabilístics assumeixen una varietat de models de dades i apliquen l’estimació de màxima versemblança i els criteris de Bayes per identificar el model més probable i el nombre de grups.

5.1 Estudi de customers

Agafem el model name VVV, (multivariate mixture, ellipsoidal, varying volume, shape, and orientation). És el nostre cas.

sample <- scale(customer.data.clean, center=TRUE, scale=TRUE)
mclust <- Mclust(sample[,1:4], modelName = "VVV")
summary(mclust)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VVV (ellipsoidal, varying volume, shape, and orientation) model with 3
## components: 
## 
##  log-likelihood    n df       BIC       ICL
##       -13584.47 3895 44 -27532.71 -28305.64
## 
## Clustering table:
##    1    2    3 
## 1340 1044 1511
plot(mclust, what = 'BIC', legendArgs = list(x = "bottomright", ncol = 10))

plot(mclust, what = "uncertainty")

Ens ha generat 3 grups.

Anem a veure amb variables 2 a 2: ### Days x TotalCost

mclust <- Mclust(sample[,c("Days","TotalCost")], modelName = "VVV")
summary(mclust)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VVV (ellipsoidal, varying volume, shape, and orientation) model with 9
## components: 
## 
##  log-likelihood    n df       BIC       ICL
##       -7188.852 3895 53 -14815.88 -17161.94
## 
## Clustering table:
##   1   2   3   4   5   6   7   8   9 
## 625 176 411 414 687 444  98 582 458
plot(mclust, what = 'BIC', legendArgs = list(x = "bottomright", ncol = 5), modelName = "VVV")

plot(mclust, what = "uncertainty")

Amb la l’estudi de Days i total cost ens està mostrant fins a 9 grups diferents. De 4 a 9 grups el BIC (Bayesian Information Criterion) és molt similar.

5.1.1 Days x NumStocks

mclust <- Mclust(sample[,c("Days","NumStocks")], modelName = "VVV")
summary(mclust)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VVV (ellipsoidal, varying volume, shape, and orientation) model with 9
## components: 
## 
##  log-likelihood    n df       BIC       ICL
##       -7123.542 3895 53 -14685.26 -17379.11
## 
## Clustering table:
##   1   2   3   4   5   6   7   8   9 
## 415 422 602 743 510 498 338  99 268
plot(mclust, what = 'BIC', legendArgs = list(x = "bottomright", ncol = 5), modelName = "VVV")

plot(mclust, what = "uncertainty")

Days x NumStocks ens passa exactament el mateix.

Anem a fer un darrer estudi del BIC

BIC <- mclustBIC(sample)
plot(BIC)

summary(BIC)
## Best BIC values:
##             VEV,9       VEV,8      VEV,7
## BIC      -20733.6 -21255.9952 -23332.892
## BIC diff      0.0   -522.3999  -2599.297
mod1 <- Mclust(sample, x = BIC)
summary(mod1, parameters = TRUE)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VEV (ellipsoidal, equal shape) model with 9 components: 
## 
##  log-likelihood    n  df      BIC       ICL
##       -9912.088 3895 110 -20733.6 -21589.19
## 
## Clustering table:
##   1   2   3   4   5   6   7   8   9 
## 863 162 329 675 152 298 710 592 114 
## 
## Mixing probabilities:
##          1          2          3          4          5          6          7 
## 0.21259164 0.04382844 0.08530769 0.16654909 0.04180480 0.08144189 0.18138048 
##          8          9 
## 0.15040885 0.03668712 
## 
## Means:
##                   [,1]        [,2]         [,3]       [,4]         [,5]
## Days        -0.8216975 -0.09640612  0.003541239 -0.1265173  0.326200610
## NumInvoices  0.8623392  0.32540992  0.105643460 -0.3588188  0.005405116
## TotalCost    0.8399817  0.52682303 -0.069819382 -0.4124888 -0.443884545
## NumStocks    0.6317835  0.03190393 -0.119282024 -0.3745369  0.440404839
##                   [,6]       [,7]       [,8]      [,9]
## Days        -0.4134963 -0.3295983  1.7172975 0.5780150
## NumInvoices  1.1834385 -0.7475808 -0.7475808 0.1251766
## TotalCost    1.3064223 -0.7232285 -0.7645318 0.8538076
## NumStocks    1.2386544 -0.5842337 -0.6478521 0.5714852
## 
## Variances:
## [,,1]
##                     Days NumInvoices   TotalCost   NumStocks
## Days         0.006662607 -0.04057397 -0.03581171 -0.02021258
## NumInvoices -0.040573972  1.68784994  1.26853213  0.42003817
## TotalCost   -0.035811706  1.26853213  1.71647894  0.60608890
## NumStocks   -0.020212585  0.42003817  0.60608890  1.34835971
## [,,2]
##                   Days NumInvoices  TotalCost  NumStocks
## Days         0.8626668  -0.2272808 -0.1330553 -0.2681446
## NumInvoices -0.2272808   0.5775007  0.4543723  0.4742551
## TotalCost   -0.1330553   0.4543723  1.1673724  1.0738690
## NumStocks   -0.2681446   0.4742551  1.0738690  1.0224531
## [,,3]
##                   Days NumInvoices  TotalCost   NumStocks
## Days         1.3401967 -0.25010132 -0.2553406 -0.24003980
## NumInvoices -0.2501013  0.24985736  0.0811263  0.05356166
## TotalCost   -0.2553406  0.08112630  0.3645501  0.33102521
## NumStocks   -0.2400398  0.05356166  0.3310252  0.30788137
## [,,4]
##                      Days   NumInvoices     TotalCost     NumStocks
## Days         4.125584e-01  4.437192e-20 -1.033005e-02 -1.940712e-02
## NumInvoices  4.437192e-20  7.125119e-04  2.645351e-17 -1.673268e-18
## TotalCost   -1.033005e-02  2.645351e-17  1.017483e-01  4.379301e-02
## NumStocks   -1.940712e-02 -1.673268e-18  4.379301e-02  9.161723e-02
## [,,5]
##                   Days NumInvoices  TotalCost  NumStocks
## Days         2.7065908  -0.3722718 -0.3036071 -0.5224246
## NumInvoices -0.3722718   0.9301467  0.4021086  0.1401498
## TotalCost   -0.3036071   0.4021086  0.2359166  0.2183954
## NumStocks   -0.5224246   0.1401498  0.2183954  0.4710598
## [,,6]
##                   Days NumInvoices  TotalCost  NumStocks
## Days         0.1367154  -0.6524757 -0.1597707 -0.9568664
## NumInvoices -0.6524757   7.3456921  2.0439347  3.8724098
## TotalCost   -0.1597707   2.0439347  2.4725438  1.9786256
## NumStocks   -0.9568664   3.8724098  1.9786256  9.4840355
## [,,7]
##                      Days   NumInvoices    TotalCost     NumStocks
## Days         1.579166e-01 -2.116084e-19 1.602121e-03 -7.442771e-03
## NumInvoices -2.116084e-19  2.724293e-04 1.507792e-17  1.182417e-17
## TotalCost    1.602121e-03  1.507792e-17 3.116479e-02  1.558987e-02
## NumStocks   -7.442771e-03  1.182417e-17 1.558987e-02  4.259392e-02
## [,,8]
##                      Days   NumInvoices     TotalCost     NumStocks
## Days         1.738108e-01 -1.054055e-19 -1.566002e-03 -2.286618e-03
## NumInvoices -1.054055e-19  2.990525e-04 -7.976302e-17  4.282695e-17
## TotalCost   -1.566002e-03 -7.976302e-17  3.644557e-02  1.784302e-02
## NumStocks   -2.286618e-03  4.282695e-17  1.784302e-02  4.405941e-02
## [,,9]
##                  Days NumInvoices  TotalCost NumStocks
## Days         2.429545  -1.2358349 -1.7067809 -1.634157
## NumInvoices -1.235835   0.6727232  0.8038324  1.139690
## TotalCost   -1.706781   0.8038324  2.7577851  1.726544
## NumStocks   -1.634157   1.1396904  1.7265439  5.955426

5.2 Estudi de productes

sample <- scale(product.data.clean, center=TRUE, scale=TRUE)
mclust <- Mclust(sample[,1:3], modelName = "VVV")
summary(mclust)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VVV (ellipsoidal, varying volume, shape, and orientation) model with 9
## components: 
## 
##  log-likelihood    n df       BIC       ICL
##       -2062.095 3380 89 -4847.372 -6161.645
## 
## Clustering table:
##   1   2   3   4   5   6   7   8   9 
## 585 350 608 450 228 173 584 167 235
plot(mclust, what = 'BIC', legendArgs = list(x = "bottomright", ncol = 5), modelName = "VVV")

plot(mclust, what = "uncertainty")

9 clusters en total. I de 6 a 9 clusters valors molt similars.

5.2.1 UnitPrice x Days

mclust <- Mclust(sample[,c("Days","UnitPrice")], modelName = "VVV")
summary(mclust)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VVV (ellipsoidal, varying volume, shape, and orientation) model with 8
## components: 
## 
##  log-likelihood    n df      BIC       ICL
##       -2072.217 3380 47 -4526.34 -5710.032
## 
## Clustering table:
##   1   2   3   4   5   6   7   8 
## 725 516 525 261 409  92 454 398
plot(mclust, what = 'BIC', legendArgs = list(x = "bottomright", ncol = 5), modelName = "VVV")

plot(mclust, what = "uncertainty")

UnitPrice x Days ens dóna 9 clusters, pero de 4 a 8 és molt similar.

5.2.2 UnitPrice x TotalStockUnits

mclust <- Mclust(sample[,c("UnitPrice","TotalStockUnits")], modelName = "VVV")
summary(mclust)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VVV (ellipsoidal, varying volume, shape, and orientation) model with 6
## components: 
## 
##  log-likelihood    n df       BIC      ICL
##        -3923.13 3380 35 -8130.657 -9556.23
## 
## Clustering table:
##   1   2   3   4   5   6 
## 950 190 524 463 536 717
plot(mclust, what = 'BIC', legendArgs = list(x = "bottomright", ncol = 5), modelName = "VVV")

plot(mclust, what = "uncertainty")

UnitPrice x TotalStockUnits de 7 a 9 clusters

5.2.3 TotalStockUnits x Days

mclust <- Mclust(sample[,c("TotalStockUnits","Days")], modelName = "VVV")
summary(mclust)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VVV (ellipsoidal, varying volume, shape, and orientation) model with 8
## components: 
## 
##  log-likelihood    n df      BIC      ICL
##        2519.101 3380 47 4656.297 3341.915
## 
## Clustering table:
##   1   2   3   4   5   6   7   8 
## 234 632 128 375 412 894 422 283
plot(mclust, what = 'BIC', legendArgs = list(x = "bottomright", ncol = 5), modelName = "VVV")

plot(mclust, what = "uncertainty")

Anem a fer un altre estudi del BIC

BIC <- mclustBIC(sample)
plot(BIC)

summary(BIC)
## Best BIC values:
##              VVE,9      VVV,9      VVI,9
## BIC      -3919.759 -4396.0839 -4483.5522
## BIC diff     0.000  -476.3252  -563.7936
mod1 <- Mclust(sample, x = BIC)
summary(mod1, parameters = TRUE)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VVE (ellipsoidal, equal orientation) model with 9 components: 
## 
##  log-likelihood    n df       BIC       ICL
##       -1695.796 3380 65 -3919.759 -5298.618
## 
## Clustering table:
##   1   2   3   4   5   6   7   8   9 
## 566 486  70 459 435 278 266 466 354 
## 
## Mixing probabilities:
##          1          2          3          4          5          6          7 
## 0.14716265 0.14459761 0.02238242 0.14674063 0.12349818 0.09292232 0.08689564 
##          8          9 
## 0.12104182 0.11475873 
## 
## Means:
##                       [,1]       [,2]       [,3]       [,4]       [,5]
## UnitPrice       -0.7315634  0.5734528 -0.7112881  0.8292693  0.2639911
## Days            -0.5036560  1.3779237  0.1362377 -0.5063829 -0.4223195
## TotalStockUnits  0.1180646 -0.4652270  2.2871462 -0.2564182 -0.4207641
##                       [,6]       [,7]        [,8]       [,9]
## UnitPrice       -0.4865466 -0.2410860 -0.28897304 -0.1088588
## Days            -0.3263965 -0.5200799 -0.51688709  1.1883656
## TotalStockUnits -0.1603981  1.7290312 -0.02533168 -0.3832369
## 
## Variances:
## [,,1]
##                     UnitPrice          Days TotalStockUnits
## UnitPrice        7.293574e-02 -5.268068e-05    8.367072e-05
## Days            -5.268068e-05  2.517675e-04    6.261126e-05
## TotalStockUnits  8.367072e-05  6.261126e-05    1.812563e-01
## [,,2]
##                     UnitPrice          Days TotalStockUnits
## UnitPrice        1.3439005606  0.0003799178   -1.037973e-03
## Days             0.0003799178  1.8682837538   -6.461269e-04
## TotalStockUnits -0.0010379728 -0.0006461269    1.305055e-05
## [,,3]
##                    UnitPrice         Days TotalStockUnits
## UnitPrice       0.1581333816 0.0002837577     0.011152469
## Days            0.0002837577 0.5440980293     0.004858508
## TotalStockUnits 0.0111524691 0.0048585084    14.599448818
## [,,4]
##                     UnitPrice          Days TotalStockUnits
## UnitPrice        1.2136916781 -8.803005e-04   -9.163063e-04
## Days            -0.0008803005  2.385458e-04    9.856040e-06
## TotalStockUnits -0.0009163063  9.856040e-06    2.678426e-02
## [,,5]
##                     UnitPrice          Days TotalStockUnits
## UnitPrice        0.8506246414 -6.132157e-04   -6.557828e-04
## Days            -0.0006132157  5.343750e-03   -9.640934e-07
## TotalStockUnits -0.0006557828 -9.640934e-07    1.185384e-03
## [,,6]
##                     UnitPrice          Days TotalStockUnits
## UnitPrice        0.2681562487 -1.782899e-04   -1.661928e-04
## Days            -0.0001782899  2.238269e-02    1.067855e-05
## TotalStockUnits -0.0001661928  1.067855e-05    5.287550e-02
## [,,7]
##                     UnitPrice          Days TotalStockUnits
## UnitPrice        0.6992368420 -5.067393e-04    0.0009117609
## Days            -0.0005067393  3.346052e-05    0.0006501426
## TotalStockUnits  0.0009117609  6.501426e-04    1.8796373845
## [,,8]
##                     UnitPrice          Days TotalStockUnits
## UnitPrice        0.3194221003 -2.316995e-04   -1.696151e-04
## Days            -0.0002316995  1.040126e-06    3.463848e-05
## TotalStockUnits -0.0001696151  3.463848e-05    9.968597e-02
## [,,9]
##                     UnitPrice          Days TotalStockUnits
## UnitPrice        0.8546688172  0.0001050217   -0.0006549083
## Days             0.0001050217  0.9998011482   -0.0003433872
## TotalStockUnits -0.0006549083 -0.0003433872    0.0066835656

Els millors BIC (Bayesian Information Criterion) que tenim és de 9 clusters.


6 Clustering: Comparació de resultats


Partim de la base que les dades que tenim són de comportament de compra dels clients i de tipologia de productes:

6.1 Customers

Fem un resum dels resultats obtinguts amb els 3 tipus de fer clustering pels clients.

6.1.1 K-Means

Hem decidit que teníem entre 2 i 3 clusters, però en algun dels estudis hem tingut fins a 7 possibles clusters.

6.1.2 Hierarchical Agglomerative Clustering (HAC)

Són els que ens han donat uns resultats més clars. Encara que amb enllaç complet i promig hem tingut de 3 a 7 clusters amb mètode ward hem pogut fixar millor el nombre de clusters a 4 o 5.

6.1.3 Mètodes d’agregació probabilístics

En aquest cas hem tingut més disparitat de dades. En el primer estudi ens ha donat 3 clusters i després 5,6 i 9. No m’ha generat molta fiabilitat l’estudi fet amb mètodes probabilistics.

De totes maneres hem vist en tots els estudis que hi havia un ventall una mica gran de possibles tipologies de clients, però si hagués de fer un estudi, seria entre 3 i 5 clusters.

6.2 Products

Fem el mateix resum però pels productes.

6.2.1 K-Means

En aquest cas hem tingut molt clarament una agrupació entre 4 i 5 clusters i més concretament en 4.

6.2.2 Hierarchical Agglomerative Clustering (HAC)

Mateixos resultats que amb K-Means, 4 i 5 clusters.

6.2.3 Mètodes d’agregació probabilístics

Un altre cop tenim valors molt dispars comparats amb els mètodes basats en les distàncies HAC i K-Means.


7 Classificació: Model supervisat sense haver aplicant prèviament PCA/SVD


Per a crear el model de l’arbre de decisió anem a fer servir la llibreria C50. Aquesta llibreria implementa funcions per a generar arbres amb l’algoritme C5.0 alhora basat en l’algoritme ID3 de Quinlan. Com a millores, permet treballar amb variables categòriques i discretes i realitza la poda de forma automàtica.

El mètode ID3 tracta de trobar una partició que asseguri la màxima capacitat predictiva i la màxima homogeneïtat de les classes. Volem pocs atributs per a predir el màxim possible. Per a fer-ho, a l’escollir un atribut, calcula el guany obtingut, i el càlcul el fa amb la diferència de la informació d’una partició respecte al desordre (o entropia) que es genera a partir d’aquella elecció.

La millora en la poda la realitza calculant la taxa d’error de les prediccions de les fulles. Es realitza una estimació pessimista de la predicció (amb calcul de la distribucuó binomial), i si els nodes fills tenen una estimació pitjor que la dels pares s’eliminen. És postpoda, és a dir, es realitza la poda després de contruir l’arbre.

En aquest primer estudi anem a executar l’algoritme de predicció sense excloure cap dels atributs del dataset.

7.1 Preparació de les dades

Prepararem les dades per a avaluar l’arbre de decisió. Dividim en conjunt de prova i conjunt d’entrenament. Un el fem servir per a construir el model de l’arbre de decisió i el de prova per avaluar quina precisió obtenim d’aquest arbre. Fem servir la proporció 0.7 entrenament, 0.3 d’avaluació. La variable a avaluar en aquest cas és la columna ‘y’.

# separem el 70% de conjunt d'entrenament del 30% per a avaluar

## obtenim primer una mostra del total de forma aleatoria (hem guardat físicament aquest random per poder reproduïr el mateix i el carreguem de disc)

# random.data <- sample(1:nrow(bank.data), 0.70 * nrow(bank.data))
# write.csv(random.data,"random-data.csv",row.names=FALSE)

bank.data <- read.csv('bank-clean-discret.csv', sep = ',')
summary(bank.data)
##       age                 job           marital          education    
##  Min.   :1.000   blue-collar:8507   divorced: 4494   primary  : 5926  
##  1st Qu.:2.000   management :7896   married :22762   secondary:20981  
##  Median :3.000   technician :6615   single  :10997   tertiary :11346  
##  Mean   :3.046   admin.     :4565                                     
##  3rd Qu.:4.000   services   :3708                                     
##  Max.   :5.000   retired    :1492                                     
##                  (Other)    :5470                                     
##  default        balance        housing      loan            contact     
##  no :37485   Min.   :-1884.0   no :16220   no :31548   cellular :24919  
##  yes:  768   1st Qu.:   43.0   yes:22033   yes: 6705   telephone: 2169  
##              Median :  338.0                           unknown  :11165  
##              Mean   :  623.6                                            
##              3rd Qu.:  957.0                                            
##              Max.   : 3388.0                                            
##                                                                         
##       day            month          duration        campaign     
##  Min.   : 1.00   may    :12161   Min.   :1.000   Min.   : 1.000  
##  1st Qu.: 8.00   jul    : 6141   1st Qu.:2.000   1st Qu.: 1.000  
##  Median :16.00   aug    : 5320   Median :3.000   Median : 2.000  
##  Mean   :15.78   jun    : 4333   Mean   :3.003   Mean   : 2.776  
##  3rd Qu.:21.00   nov    : 2940   3rd Qu.:4.000   3rd Qu.: 3.000  
##  Max.   :31.00   apr    : 2450   Max.   :5.000   Max.   :58.000  
##                  (Other): 4908                                   
##      pdays          previous           poutcome       y        
##  Min.   : -1.0   Min.   :  0.0000   failure: 4099   no :34124  
##  1st Qu.: -1.0   1st Qu.:  0.0000   other  : 1551   yes: 4129  
##  Median : -1.0   Median :  0.0000   success: 1156              
##  Mean   : 40.3   Mean   :  0.5659   unknown:31447              
##  3rd Qu.: -1.0   3rd Qu.:  0.0000                              
##  Max.   :871.0   Max.   :275.0000                              
## 
random.data <- read.csv('random-data.csv', sep = ',')
random.data <- unlist(random.data, use.names=FALSE)

data.train <- bank.data[random.data,]
data.test <- bank.data[-random.data,]

# en les dades d'entrenament, separem les dades per a generar l'arbre de la dada objectiu (target)
train.x <- data.train[,1:16]
train.y <- data.train[,17]

# en les dades de test, separem les dades per a generar l'arbre de la dada objectiu (target)
test.x <- data.test[,1:16]
test.y <- data.test[,17]

7.2 Creació del model del arbre de decissió

model1 <- C50::C5.0(train.x, train.y, rules=TRUE)
predicted.model1 <- predict(model1, test.x, type="class")
cat(sprintf("La precissió de l'arbre és del: %.2f %%",100*sum(predicted.model1 == test.y) / length(predicted.model1)))
## La precissió de l'arbre és del: 91.20 %
summary(model1)
## 
## Call:
## C5.0.default(x = train.x, y = train.y, rules = TRUE)
## 
## 
## C5.0 [Release 2.07 GPL Edition]      Tue Jun 16 22:30:08 2020
## -------------------------------
## 
## Class specified by attribute `outcome'
## 
## Read 26777 cases (17 attributes) from undefined.data
## 
## Rules:
## 
## Rule 1: (3852/22, lift 1.1)
##  contact in {telephone, unknown}
##  duration <= 2
##  ->  class no  [0.994]
## 
## Rule 2: (152, lift 1.1)
##  housing = yes
##  month = jan
##  pdays > 155
##  ->  class no  [0.994]
## 
## Rule 3: (1168/10, lift 1.1)
##  housing = no
##  contact = cellular
##  day > 7
##  month = aug
##  duration <= 2
##  ->  class no  [0.991]
## 
## Rule 4: (7043/79, lift 1.1)
##  month in {apr, dec, jan, jul, jun, may}
##  duration <= 2
##  ->  class no  [0.989]
## 
## Rule 5: (782/11, lift 1.1)
##  contact = cellular
##  month = nov
##  duration <= 2
##  ->  class no  [0.985]
## 
## Rule 6: (9851/159, lift 1.1)
##  duration <= 2
##  previous <= 2
##  poutcome in {failure, other, unknown}
##  ->  class no  [0.984]
## 
## Rule 7: (9138/155, lift 1.1)
##  housing = yes
##  duration <= 3
##  ->  class no  [0.983]
## 
## Rule 8: (4770/112, lift 1.1)
##  day > 20
##  duration <= 3
##  pdays <= 385
##  poutcome in {failure, other, unknown}
##  ->  class no  [0.976]
## 
## Rule 9: (14615/367, lift 1.1)
##  balance <= 2174
##  duration <= 3
##  poutcome in {failure, other, unknown}
##  ->  class no  [0.975]
## 
## Rule 10: (7676/291, lift 1.1)
##  contact = unknown
##  month in {jan, jul, jun, may}
##  poutcome in {failure, other, unknown}
##  ->  class no  [0.962]
## 
## Rule 11: (930/45, lift 1.1)
##  housing = yes
##  day > 17
##  month = nov
##  ->  class no  [0.951]
## 
## Rule 12: (2896/141, lift 1.1)
##  housing = yes
##  campaign > 3
##  previous <= 5
##  ->  class no  [0.951]
## 
## Rule 13: (356/20, lift 1.1)
##  housing = yes
##  day <= 4
##  month = feb
##  previous <= 3
##  ->  class no  [0.941]
## 
## Rule 14: (3947/241, lift 1.1)
##  education = primary
##  pdays <= 385
##  previous <= 4
##  poutcome in {failure, other, unknown}
##  ->  class no  [0.939]
## 
## Rule 15: (15016/970, lift 1.0)
##  housing = yes
##  pdays <= 385
##  poutcome in {failure, other, unknown}
##  ->  class no  [0.935]
## 
## Rule 16: (18970/1323, lift 1.0)
##  month in {aug, jan, jul, may, nov}
##  pdays <= 385
##  poutcome in {failure, other, unknown}
##  ->  class no  [0.930]
## 
## Rule 17: (1062/73, lift 1.0)
##  marital in {divorced, married}
##  day <= 9
##  month in {apr, feb}
##  pdays <= 385
##  poutcome in {failure, other, unknown}
##  ->  class no  [0.930]
## 
## Rule 18: (708/49, lift 1.0)
##  day <= 20
##  month in {apr, feb}
##  pdays <= 385
##  poutcome = failure
##  ->  class no  [0.930]
## 
## Rule 19: (1357/99, lift 1.0)
##  day <= 9
##  month = feb
##  pdays <= 385
##  poutcome in {failure, other, unknown}
##  ->  class no  [0.926]
## 
## Rule 20: (15, lift 8.7)
##  housing = no
##  month = feb
##  pdays <= 92
##  poutcome = success
##  ->  class yes  [0.941]
## 
## Rule 21: (24/2, lift 8.1)
##  age > 3
##  housing = no
##  contact = cellular
##  month = feb
##  poutcome = success
##  ->  class yes  [0.885]
## 
## Rule 22: (23/2, lift 8.1)
##  marital = divorced
##  day > 17
##  duration > 2
##  campaign <= 3
##  poutcome = success
##  ->  class yes  [0.880]
## 
## Rule 23: (37/5, lift 7.8)
##  day <= 7
##  month = aug
##  poutcome = success
##  ->  class yes  [0.846]
## 
## Rule 24: (10/1, lift 7.7)
##  marital = single
##  housing = no
##  day <= 9
##  month = apr
##  duration > 3
##  ->  class yes  [0.833]
## 
## Rule 25: (93/17, lift 7.5)
##  duration > 3
##  pdays > 385
##  ->  class yes  [0.811]
## 
## Rule 26: (7/1, lift 7.2)
##  duration > 2
##  campaign > 3
##  previous > 5
##  poutcome = success
##  ->  class yes  [0.778]
## 
## Rule 27: (47/10, lift 7.1)
##  housing = no
##  day > 9
##  day <= 20
##  month = feb
##  duration > 1
##  poutcome in {other, unknown}
##  ->  class yes  [0.776]
## 
## Rule 28: (20/4, lift 7.1)
##  month = jun
##  duration > 2
##  pdays <= 385
##  previous > 4
##  poutcome in {failure, other}
##  ->  class yes  [0.773]
## 
## Rule 29: (434/99, lift 7.1)
##  housing = no
##  duration > 2
##  poutcome = success
##  ->  class yes  [0.771]
## 
## Rule 30: (429/98, lift 7.1)
##  duration > 3
##  campaign <= 3
##  poutcome = success
##  ->  class yes  [0.770]
## 
## Rule 31: (374/86, lift 7.1)
##  month in {apr, aug, dec, jul, jun, mar, oct, sep}
##  duration > 2
##  campaign <= 3
##  poutcome = success
##  ->  class yes  [0.769]
## 
## Rule 32: (154/36, lift 7.0)
##  contact = cellular
##  month in {mar, oct, sep}
##  poutcome = success
##  ->  class yes  [0.763]
## 
## Rule 33: (187/50, lift 6.7)
##  day > 20
##  month in {apr, feb}
##  duration > 3
##  ->  class yes  [0.730]
## 
## Rule 34: (5/1, lift 6.6)
##  education = primary
##  month in {mar, oct, sep}
##  duration > 1
##  previous > 4
##  ->  class yes  [0.714]
## 
## Rule 35: (195/63, lift 6.2)
##  contact in {cellular, telephone}
##  month = jun
##  duration > 3
##  pdays <= 385
##  ->  class yes  [0.675]
## 
## Rule 36: (608/235, lift 5.6)
##  education in {secondary, tertiary}
##  month in {dec, mar, oct, sep}
##  duration > 2
##  ->  class yes  [0.613]
## 
## Rule 37: (289/120, lift 5.4)
##  contact in {cellular, telephone}
##  month = jun
##  duration > 2
##  pdays <= 385
##  ->  class yes  [0.584]
## 
## Rule 38: (204/93, lift 5.0)
##  housing = no
##  day > 9
##  day <= 20
##  month in {apr, feb}
##  duration > 1
##  ->  class yes  [0.544]
## 
## Default class: no
## 
## 
## Evaluation on training data (26777 cases):
## 
##          Rules     
##    ----------------
##      No      Errors
## 
##      38 2360( 8.8%)   <<
## 
## 
##     (a)   (b)    <-classified as
##    ----  ----
##   23491   374    (a): class no
##    1986   926    (b): class yes
## 
## 
##  Attribute usage:
## 
##   97.21% poutcome
##   95.47% month
##   90.76% pdays
##   63.94% duration
##   63.44% housing
##   54.58% balance
##   51.17% previous
##   40.88% contact
##   31.44% day
##   17.03% education
##   12.87% campaign
##    4.09% marital
##    0.09% age
## 
## 
## Time: 0.8 secs

L’arbre de decissió ha fet servir un total de 13 atributs per a generar les regles i ha generat un total de 38 regles.

7.3 Visualització de l’arbre de decissió

set.seed(123)
model1 <- C50::C5.0(train.x, train.y, tree=TRUE)
summary(model1)
## 
## Call:
## C5.0.default(x = train.x, y = train.y, tree = TRUE)
## 
## 
## C5.0 [Release 2.07 GPL Edition]      Tue Jun 16 22:30:13 2020
## -------------------------------
## 
## Class specified by attribute `outcome'
## 
## Read 26777 cases (17 attributes) from undefined.data
## 
## Decision tree:
## 
## poutcome = success:
## :...duration <= 2:
## :   :...contact in {telephone,unknown}: no (13/1)
## :   :   contact = cellular:
## :   :   :...housing = yes: no (55/8)
## :   :       housing = no:
## :   :       :...month in {apr,dec,jan,jul,jun,may}: no (40/7)
## :   :           month in {mar,oct,sep}: yes (17/6)
## :   :           month = aug:
## :   :           :...day <= 7: yes (7/1)
## :   :           :   day > 7: no (19/2)
## :   :           month = nov:
## :   :           :...balance <= 719: yes (4)
## :   :           :   balance > 719: no (7/1)
## :   :           month = feb:
## :   :           :...pdays <= 92: yes (6)
## :   :               pdays > 92:
## :   :               :...age <= 3: no (5)
## :   :                   age > 3: yes (4/1)
## :   duration > 2:
## :   :...housing = no: yes (434/99)
## :       housing = yes:
## :       :...campaign > 3:
## :           :...previous <= 5: no (11/1)
## :           :   previous > 5: yes (5/1)
## :           campaign <= 3:
## :           :...month in {apr,aug,dec,jul,jun,mar,oct,sep}: yes (97/22)
## :               month = jan:
## :               :...pdays <= 155: yes (2)
## :               :   pdays > 155: no (3)
## :               month = may:
## :               :...duration <= 3: no (8/2)
## :               :   duration > 3: yes (52/22)
## :               month = nov:
## :               :...day <= 17: yes (10)
## :               :   day > 17:
## :               :   :...marital = divorced: yes (2)
## :               :       marital in {married,single}: no (10)
## :               month = feb:
## :               :...duration <= 3: no (4)
## :                   duration > 3:
## :                   :...previous > 3: yes (7)
## :                       previous <= 3:
## :                       :...day <= 4: no (4)
## :                           day > 4: yes (2)
## poutcome in {failure,other,unknown}:
## :...pdays > 385:
##     :...duration <= 3: no (58/12)
##     :   duration > 3:
##     :   :...previous <= 4: yes (63/10)
##     :       previous > 4: no (10/3)
##     pdays <= 385:
##     :...month in {aug,jan,jul,jun,may,nov}:
##         :...duration <= 2: no (8999/54)
##         :   duration > 2:
##         :   :...contact = unknown:
##         :       :...month in {jan,jul,jun,may}: no (4676/286)
##         :       :   month in {aug,nov}:
##         :       :   :...day <= 15: yes (7)
##         :       :       day > 15: no (14/1)
##         :       contact in {cellular,telephone}:
##         :       :...month in {aug,jan,jul,may,nov}: no (7988/1109)
##         :           month = jun:
##         :           :...previous > 4: yes (20/4)
##         :               previous <= 4:
##         :               :...duration <= 3: no (75/25)
##         :                   duration > 3: yes (143/52)
##         month in {dec,mar,oct,sep}:
##         :...duration <= 1: no (123/5)
##         :   duration > 1:
##         :   :...education = primary:
##         :       :...previous <= 4: no (53/16)
##         :       :   previous > 4: yes (4/1)
##         :       education in {secondary,tertiary}:
##         :       :...duration > 2:
##         :           :...poutcome = failure: no (84/35)
##         :           :   poutcome in {other,unknown}: yes (335/139)
##         :           duration <= 2:
##         :           :...previous > 2: no (18/1)
##         :               previous <= 2:
##         :               :...marital = divorced: yes (10/2)
##         :                   marital in {married,single}: no (120/40)
##         month in {apr,feb}:
##         :...day > 20:
##             :...duration <= 3: no (171/29)
##             :   duration > 3: yes (152/47)
##             day <= 20:
##             :...housing = yes: no (1849/125)
##                 housing = no:
##                 :...day > 9:
##                     :...poutcome = failure: no (25/1)
##                     :   poutcome in {other,unknown}:
##                     :   :...duration <= 1: no (17/2)
##                     :       duration > 1:
##                     :       :...month = apr: no (101/43)
##                     :           month = feb: yes (47/10)
##                     day <= 9:
##                     :...month = feb: no (707/64)
##                         month = apr:
##                         :...loan = yes: yes (6/1)
##                             loan = no:
##                             :...marital in {divorced,married}: no (44/8)
##                                 marital = single:
##                                 :...duration > 3: yes (9/1)
##                                     duration <= 3:
##                                     :...balance <= 2174: no (17/2)
##                                         balance > 2174: yes (4/1)
## 
## 
## Evaluation on training data (26777 cases):
## 
##      Decision Tree   
##    ----------------  
##    Size      Errors  
## 
##      58 2303( 8.6%)   <<
## 
## 
##     (a)   (b)    <-classified as
##    ----  ----
##   23445   420    (a): class no
##    1883  1029    (b): class yes
## 
## 
##  Attribute usage:
## 
##  100.00% poutcome
##   97.58% month
##   96.98% pdays
##   90.17% duration
##   48.92% contact
##   13.60% housing
##   12.04% day
##    2.33% education
##    2.04% previous
##    0.81% campaign
##    0.81% marital
##    0.30% loan
##    0.12% balance
##    0.03% age
## 
## 
## Time: 0.5 secs

L’arbre té un tamany de 58 nodes amb una taxa d’error del 8.6%

7.4 Validació de la qualitat del model

predicted.model1 <- predict(model1, test.x, type="class")
print(sprintf("La precissió de l'arbre és del: %.2f %%",100*sum(predicted.model1 == test.y) / length(predicted.model1)))
## [1] "La precissió de l'arbre és del: 91.24 %"

Tenim una precisió de model força bona, de més del 91% dels casos. Si per exemple, tinguéssim uns comercials que fessin campanyes telefòniques per a captar nous clients, podríem fer servir aquest model per a predir amb més d’un 90% de fiabilitat quin target atacar o quines decisions de comportament prendre.

Fent servir la matriu de confusió veiem la qualitat del model també:

## creem la matriu per avaluar la prova
matrix.conf <- table(Class=predicted.model1,Predicted=test.y)
percent.correct <- 100 * sum(diag(matrix.conf)) / sum(matrix.conf)
print(sprintf("L'error de classificació és: %.2f %%",100 - percent.correct))
## [1] "L'error de classificació és: 8.76 %"
mosaicplot(matrix.conf)

confusionMatrix(predicted.model1,test.y, dnn = c("Prediction"))
## Confusion Matrix and Statistics
## 
##           NA
## Prediction    no   yes
##        no  10066   812
##        yes   193   405
##                                           
##                Accuracy : 0.9124          
##                  95% CI : (0.9071, 0.9175)
##     No Information Rate : 0.894           
##     P-Value [Acc > NIR] : 2.267e-11       
##                                           
##                   Kappa : 0.4047          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.9812          
##             Specificity : 0.3328          
##          Pos Pred Value : 0.9254          
##          Neg Pred Value : 0.6773          
##              Prevalence : 0.8940          
##          Detection Rate : 0.8771          
##    Detection Prevalence : 0.9479          
##       Balanced Accuracy : 0.6570          
##                                           
##        'Positive' Class : no              
## 

La veritat és que sense fer PCA tenim un valor de predicció molt bó amb un Kappa (aletorietat de la predicció) realment alt, cosa que es bó.

Fleiss et al. (2003) va afirmar que, a la majoria dels propòsits,

  • És poden prendre valors superiors a 0,75 per representar un excel·lent acord fora de l’atzar,
  • És poden prendre valors inferiors a 0,40 per representar un mal acord fora de l’atzar i
  • És poden prendre valors entre 0,40 i 0,75 per representar un acord just per a un bon acord fora de l’atzar.

Altre cosa molt interessant:

  • Predicció del ‘no’: 0.9254
  • Predicció del ‘yes’: 0.6773

Predim molt millor el ‘no’ que el ‘yes’. Aquesta part no m’agrada tant.


8 Classificació: Model supervisat aplicant prèviament PCA/SVD


8.1 Preparació de les dades

Preparem altre cop les dades per a avaluar l’arbre de decisió. Dividim en conjunt de prova i conjunt d’entrenament. Un el fem servir per a construir el model de l’arbre de decisió i el de prova per avaluar quina precisió obtenim d’aquest arbre. Fem servir la proporció 0.7 entrenament, 0.3 d’avaluació. La variable a avaluar en aquest cas és la columna ‘y’.

bank.data <- read.csv('bank-clean-discret.csv', sep = ',')
summary(bank.data)
##       age                 job           marital          education    
##  Min.   :1.000   blue-collar:8507   divorced: 4494   primary  : 5926  
##  1st Qu.:2.000   management :7896   married :22762   secondary:20981  
##  Median :3.000   technician :6615   single  :10997   tertiary :11346  
##  Mean   :3.046   admin.     :4565                                     
##  3rd Qu.:4.000   services   :3708                                     
##  Max.   :5.000   retired    :1492                                     
##                  (Other)    :5470                                     
##  default        balance        housing      loan            contact     
##  no :37485   Min.   :-1884.0   no :16220   no :31548   cellular :24919  
##  yes:  768   1st Qu.:   43.0   yes:22033   yes: 6705   telephone: 2169  
##              Median :  338.0                           unknown  :11165  
##              Mean   :  623.6                                            
##              3rd Qu.:  957.0                                            
##              Max.   : 3388.0                                            
##                                                                         
##       day            month          duration        campaign     
##  Min.   : 1.00   may    :12161   Min.   :1.000   Min.   : 1.000  
##  1st Qu.: 8.00   jul    : 6141   1st Qu.:2.000   1st Qu.: 1.000  
##  Median :16.00   aug    : 5320   Median :3.000   Median : 2.000  
##  Mean   :15.78   jun    : 4333   Mean   :3.003   Mean   : 2.776  
##  3rd Qu.:21.00   nov    : 2940   3rd Qu.:4.000   3rd Qu.: 3.000  
##  Max.   :31.00   apr    : 2450   Max.   :5.000   Max.   :58.000  
##                  (Other): 4908                                   
##      pdays          previous           poutcome       y        
##  Min.   : -1.0   Min.   :  0.0000   failure: 4099   no :34124  
##  1st Qu.: -1.0   1st Qu.:  0.0000   other  : 1551   yes: 4129  
##  Median : -1.0   Median :  0.0000   success: 1156              
##  Mean   : 40.3   Mean   :  0.5659   unknown:31447              
##  3rd Qu.: -1.0   3rd Qu.:  0.0000                              
##  Max.   :871.0   Max.   :275.0000                              
## 
# separem el 70% de conjunt d'entrenament del 30% per a avaluar

## obetnim primer una mostra del total de forma aleatoria (hem guardat físicament aquest random per poder reproduïr el mateix i el carreguem de disc)

random.data <- read.csv('random-data.csv', sep = ',')
random.data <- unlist(random.data, use.names=FALSE)


data.train <- bank.data[random.data,]
data.test <- bank.data[-random.data,]

# en les dades d'entrenament, separem les dades per a generar l'arbre de la dada objectiu (target)
train.x <- data.train[,1:16]
train.y <- data.train[,17]

# en les dades de test, separem les dades per a generar l'arbre de la dada objectiu (target)
test.x <- data.test[,1:16]
test.y <- data.test[,17]

8.2 Creació del model del ’arbre de decisió

L’executem de forma que ens generi regles per a fer prediccions. Ara executem l’algoritme amb l’opció Winnowed. Aquesta opció no farà servir aquelles dimensions que no trobi que siguin rellevants, es a dir, està fent PCA de forma transparent i automàtica.

# amb l'atribut winnow eliminem dimensions que no tenen impacte en les classes de sortida
model2 <- C50::C5.0(train.x, train.y, rules=TRUE, size=4, control = C5.0Control(winnow = TRUE))
summary(model2)
## 
## Call:
## C5.0.default(x = train.x, y = train.y, rules = TRUE, control
##  = C5.0Control(winnow = TRUE), size = 4)
## 
## 
## C5.0 [Release 2.07 GPL Edition]      Tue Jun 16 22:30:16 2020
## -------------------------------
## 
## Class specified by attribute `outcome'
## 
## Read 26777 cases (17 attributes) from undefined.data
## 
## 7 attributes winnowed
## Estimated importance of remaining attributes:
## 
##      21%  contact
##       6%  poutcome
##       3%  month
##       2%  duration
##       1%  housing
##       1%  day
##       1%  loan
##       1%  pdays
##       1%  balance
## 
## Rules:
## 
## Rule 1: (152, lift 1.1)
##  housing = yes
##  month = jan
##  pdays > 155
##  ->  class no  [0.994]
## 
## Rule 2: (4295/27, lift 1.1)
##  housing = yes
##  day <= 19
##  month = may
##  duration <= 4
##  ->  class no  [0.993]
## 
## Rule 3: (10681/228, lift 1.1)
##  duration <= 2
##  ->  class no  [0.979]
## 
## Rule 4: (7737/302, lift 1.1)
##  contact = unknown
##  month in {aug, jan, jul, jun, may, nov}
##  poutcome in {failure, other, unknown}
##  ->  class no  [0.961]
## 
## Rule 5: (930/45, lift 1.1)
##  housing = yes
##  day > 17
##  month = nov
##  ->  class no  [0.951]
## 
## Rule 6: (256/13, lift 1.1)
##  housing = yes
##  loan = yes
##  month = jun
##  ->  class no  [0.946]
## 
## Rule 7: (606/35, lift 1.1)
##  housing = yes
##  day <= 7
##  month = feb
##  ->  class no  [0.941]
## 
## Rule 8: (1183/79, lift 1.0)
##  housing = yes
##  day <= 20
##  month = apr
##  ->  class no  [0.932]
## 
## Rule 9: (25818/2325, lift 1.0)
##  pdays <= 385
##  poutcome in {failure, other, unknown}
##  ->  class no  [0.910]
## 
## Rule 10: (9/1, lift 7.5)
##  contact = cellular
##  day <= 5
##  month = dec
##  duration > 2
##  poutcome in {failure, other, unknown}
##  ->  class yes  [0.818]
## 
## Rule 11: (93/17, lift 7.5)
##  duration > 3
##  pdays > 385
##  ->  class yes  [0.811]
## 
## Rule 12: (47/10, lift 7.1)
##  housing = no
##  day > 9
##  day <= 20
##  month = feb
##  duration > 1
##  poutcome in {other, unknown}
##  ->  class yes  [0.776]
## 
## Rule 13: (13/3, lift 6.7)
##  balance > 326
##  balance <= 1059
##  contact = cellular
##  month = dec
##  duration > 2
##  poutcome in {other, unknown}
##  ->  class yes  [0.733]
## 
## Rule 14: (135/37, lift 6.6)
##  month = sep
##  duration > 3
##  ->  class yes  [0.723]
## 
## Rule 15: (88/24, lift 6.6)
##  loan = no
##  month = mar
##  duration > 2
##  poutcome in {other, unknown}
##  ->  class yes  [0.722]
## 
## Rule 16: (5/1, lift 6.6)
##  day > 3
##  day <= 4
##  month in {jun, nov}
##  duration > 2
##  poutcome = failure
##  ->  class yes  [0.714]
## 
## Rule 17: (83/24, lift 6.5)
##  contact in {cellular, telephone}
##  day <= 3
##  month in {jul, jun, nov}
##  duration > 3
##  poutcome = unknown
##  ->  class yes  [0.706]
## 
## Rule 18: (13/4, lift 6.1)
##  balance > 1784
##  month = dec
##  duration > 1
##  poutcome in {failure, other, unknown}
##  ->  class yes  [0.667]
## 
## Rule 19: (141/48, lift 6.0)
##  day > 15
##  month = oct
##  duration > 2
##  ->  class yes  [0.657]
## 
## Rule 20: (828/309, lift 5.8)
##  poutcome = success
##  ->  class yes  [0.627]
## 
## Rule 21: (657/260, lift 5.6)
##  month in {dec, mar, oct, sep}
##  duration > 2
##  ->  class yes  [0.604]
## 
## Rule 22: (161/68, lift 5.3)
##  pdays > 385
##  ->  class yes  [0.577]
## 
## Rule 23: (204/93, lift 5.0)
##  housing = no
##  day > 9
##  day <= 20
##  month in {apr, feb}
##  duration > 1
##  ->  class yes  [0.544]
## 
## Rule 24: (230/105, lift 5.0)
##  contact in {cellular, telephone}
##  day <= 4
##  month in {jul, jun, nov}
##  duration > 2
##  ->  class yes  [0.543]
## 
## Default class: no
## 
## 
## Evaluation on training data (26777 cases):
## 
##          Rules     
##    ----------------
##      No      Errors
## 
##      24 2385( 8.9%)   <<
## 
## 
##     (a)   (b)    <-classified as
##    ----  ----
##   23573   292    (a): class no
##    2093   819    (b): class yes
## 
## 
##  Attribute usage:
## 
##   99.55% poutcome
##   97.03% pdays
##   52.73% duration
##   52.62% month
##   29.83% contact
##   28.48% housing
##   28.38% day
##    1.28% loan
##    0.10% balance
## 
## 
## Time: 0.7 secs
predicted.model2 <- predict(model2, test.x, type="class")
cat(sprintf("La precissió de l'arbre és del: %.2f %%",100*sum(predicted.model2 == test.y) / length(predicted.model2)))
## La precissió de l'arbre és del: 91.06 %

Anem a veure els resultats:

  • S’han llegit 26777 casos amb 9 atributs. Amb el winnowed s’han descartat 7 atributs per a construir les regles.

  • Classe per defecte -> ‘no’ (no contracta crèdit bancari)

  • S’han generat 24 regles de classificació

  • Ha fet servir els atributs:

    • 99.55% poutcome
    • 97.03% pdays
    • 52.73% duration
    • 52.62% month
    • 29.83% contact
    • 28.48% housing
    • 28.38% day
    • 1.28% loan
    • 0.10% balance
  • L’atribut amb més pes i que és l’arrel de l’arbre es poutcome.

  • L’arbre ha classificat malament 24 de 2385 cassos (8.9% d’errors):

  • Tots els lift, aleatorietat de la predicció estan per sobre d’1, per tant les regles són prou vàlides.

8.3 Regles de decisió generades

L’arbre ha generat 24 regles de decisió. Podem observar que la fiabilitat de les primeres regles son totes les que estan orientades a la resposta ‘no’. Ens es molt més fàcil predir el ‘no’ que el ‘yes’. La primera regla amb ‘yes’, és la 10 amb un confidence de 0.818.

  • Regla 1: housing = yes & month = jan & pdays > 155 => no contracta amb un 99.4% de validesa.
    • No contracten usuaris amb prèstec habitatge, durant el mes de gener i fa més de 155 dies del contacte d’una campanya anterior.
  • Regla 2: housing = yes & day <= 19 & month = may & duration <= 4 => no contracta amb un 99.3% de validesa.
    • No contracten diposit a termini: tenen prèstec habitatge, contacte durant els primers dies del mes, mes de maig i la trucada una durada inferior o igual a 4 minuts.
  • Regla 3: duration <= 2 => no contracta amb un 97.9% de validesa.
    • No contracten diposit a termini: aquesta regla és força interessant, quan les trucades duran poc no es contracta.
  • Regla 4: contact = unknown & month in {aug, jan, jul, jun, may, nov} & poutcome in {failure, other, unknown} => no contracta amb un 0.961 de validesa.
    • No contracten diposit a termini: el contacte no és ni phone ni cellurar, és unknown no sabem com ha estat & mesos = aug, jan, jul, jun, may, nov & resultat de la campanya no és success.
  • Regla 5: housing = yes & day > 17 & month = nov => no contracta amb un 0.951 de validesa
    • No contracten diposit a termini: té prestec d’habitatge, dia de la trucada és > 17, més de novembre i resultat de la campanya no és success
  • Regla 6: housing = yes & loan = yes & month = jun => no contracta amb un 0.946 de validesa
    • No contracten diposit a termini: té prestec d’habitatge, té prestec personal i és el mes de juny
  • Regla 7: housing = yes & day <= 7 & month = feb => no contracta amb un 0.941 de validesa
    • No contracten diposit a termini: té prestec d’habitatge, dia del mes menor de 8 i és el mes de febrer
  • Regla 8: housing = yes & day <= 20 & month = apr => no contracta amb un 0.932 de validesa
    • No contracten diposit a termini: té prestec d’habitatge, dia del mes menor de 21 i és el mes d’abril
  • Regla 9: pdays <= 385 & poutcome in {failure, other, unknown} => no contracta amb un 0.910 de validesa
    • No contracten diposit a termini: han passat 385 dies del darrer contacte de la campanya anterior i no el resultat no va ser success
  • Regla 10: contact = cellular & day <= 5 & month = dec & duration > 2 & poutcome in {failure, other, unknown} => sí contracta amb un 0.818 de validesa
    • Sí contracten diposit a termini: contacte per telèfon mòbil, dia del mes inferior a 6, el mes de desembre, durada de la trucada > 2 minuts i resultat de la darrera campanya no va ser success
  • Regla 11: duration > 3 & pdays > 385 => sí contracta amb un 0.811 de validesa
    • Sí contracten diposit a termini:durada de la trucada > 3 minuts i fa més de 385 díes del darrer contacte en la darrera campanya
  • Regla 14: month = sep & duration > 3 => sí contracta amb un 0.723 de validesa
    • Sí contracten diposit a termini: mes de setembre i durada de la trucada per sobre de 3 minuts
  • Regla 15: loan = no & month = mar & duration > 2 & poutcome in {other, unknown} => sí contracta amb un 0.722 de validesa
    • Sí contracten diposit a termini: no té prestec personal, mes de març, durada de la trucada per sobre de dos minuts i resultat de la darrera campanya és ‘unknown’ o ‘other’

8.4 Visualització de l’arbre de decissió

set.seed(123)
model2 <- C50::C5.0(train.x, train.y, control = C5.0Control(winnow = TRUE))
plot(model2)

summary(model2)
## 
## Call:
## C5.0.default(x = train.x, y = train.y, control = C5.0Control(winnow = TRUE))
## 
## 
## C5.0 [Release 2.07 GPL Edition]      Tue Jun 16 22:30:19 2020
## -------------------------------
## 
## Class specified by attribute `outcome'
## 
## Read 26777 cases (17 attributes) from undefined.data
## 
## 7 attributes winnowed
## Estimated importance of remaining attributes:
## 
##      21%  contact
##       6%  poutcome
##       3%  month
##       2%  duration
##       1%  housing
##       1%  day
##       1%  loan
##       1%  pdays
##       1%  balance
## 
## Decision tree:
## 
## poutcome = success:
## :...duration <= 2: no (177/49)
## :   duration > 2:
## :   :...housing = no: yes (434/99)
## :       housing = yes:
## :       :...month in {aug,dec,jul,mar,oct,sep}: yes (67/11)
## :           month = apr:
## :           :...day <= 20: no (14/4)
## :           :   day > 20: yes (12/2)
## :           month = jan:
## :           :...pdays <= 155: yes (2)
## :           :   pdays > 155: no (3)
## :           month = jun:
## :           :...loan = no: yes (9/1)
## :           :   loan = yes: no (2)
## :           month = nov:
## :           :...day <= 17: yes (10)
## :           :   day > 17: no (12/2)
## :           month = feb:
## :           :...duration <= 3: no (4)
## :           :   duration > 3:
## :           :   :...day <= 7: no (8/3)
## :           :       day > 7: yes (6)
## :           month = may:
## :           :...duration > 4: yes (34/12)
## :               duration <= 4:
## :               :...day <= 19: no (30/8)
## :                   day > 19: yes (4)
## poutcome in {failure,other,unknown}:
## :...pdays > 385:
##     :...duration <= 3: no (58/12)
##     :   duration > 3: yes (73/17)
##     pdays <= 385:
##     :...month in {apr,feb}:
##         :...day > 20:
##         :   :...duration <= 3: no (171/29)
##         :   :   duration > 3: yes (152/47)
##         :   day <= 20:
##         :   :...housing = yes: no (1849/125)
##         :       housing = no:
##         :       :...day <= 9:
##         :           :...month = feb: no (707/64)
##         :           :   month = apr:
##         :           :   :...loan = no: no (74/21)
##         :           :       loan = yes: yes (6/1)
##         :           day > 9:
##         :           :...poutcome = failure: no (25/1)
##         :               poutcome in {other,unknown}:
##         :               :...duration <= 1: no (17/2)
##         :                   duration > 1:
##         :                   :...month = apr: no (101/43)
##         :                       month = feb: yes (47/10)
##         month in {aug,jan,jul,jun,may,nov}:
##         :...day > 4: no (20606/1395)
##         :   day <= 4:
##         :   :...contact = unknown: no (775/33)
##         :       contact in {cellular,telephone}:
##         :       :...duration <= 2: no (183/12)
##         :           duration > 2:
##         :           :...month in {aug,jan,may}: no (155/44)
##         :               month in {jul,jun,nov}:
##         :               :...duration <= 3:
##         :                   :...day <= 2: no (42/11)
##         :                   :   day > 2: yes (29/13)
##         :                   duration > 3:
##         :                   :...poutcome = other: no (7/2)
##         :                       poutcome = failure:
##         :                       :...day <= 3: no (24/10)
##         :                       :   day > 3: yes (4/1)
##         :                       poutcome = unknown:
##         :                       :...day <= 3: yes (83/24)
##         :                           day > 3: no (14/4)
##         month in {dec,mar,oct,sep}:
##         :...duration <= 1: no (123/5)
##             duration > 1:
##             :...duration <= 2: no (163/50)
##                 duration > 2:
##                 :...poutcome = failure: no (93/39)
##                     poutcome in {other,unknown}:
##                     :...month = mar:
##                         :...loan = no: yes (87/24)
##                         :   loan = yes: no (6/2)
##                         month = oct:
##                         :...day <= 15: no (44/11)
##                         :   day > 15: yes (85/30)
##                         month = sep:
##                         :...duration <= 3: no (24/8)
##                         :   duration > 3: yes (71/25)
##                         month = dec:
##                         :...contact in {telephone,unknown}: no (4)
##                             contact = cellular:
##                             :...day <= 5: yes (7/1)
##                                 day > 5:
##                                 :...balance > 1784: yes (8/3)
##                                     balance <= 1784:
##                                     :...balance > 1059: no (8/1)
##                                         balance <= 1059:
##                                         :...balance <= 326: no (13/5)
##                                             balance > 326: yes (11/3)
## 
## 
## Evaluation on training data (26777 cases):
## 
##      Decision Tree   
##    ----------------  
##    Size      Errors  
## 
##      55 2319( 8.7%)   <<
## 
## 
##     (a)   (b)    <-classified as
##    ----  ----
##   23541   324    (a): class no
##    1995   917    (b): class yes
## 
## 
##  Attribute usage:
## 
##  100.00% poutcome
##   97.23% month
##   96.93% pdays
##   94.64% day
##   12.99% housing
##   10.21% duration
##    5.11% contact
##    0.69% loan
##    0.15% balance
## 
## 
## Time: 0.5 secs

He generat un arbre de 55 nodes. En l’arbre podem veure que el primer nivell es realitza amb l’atribut poutcome. És el que aporta particions més homogènies.

8.5 Validació de la qualitat del model

predicted.model2 <- predict(model2, test.x, type="class")
print(sprintf("La precissió de l'arbre és del: %.2f %%",100*sum(predicted.model2 == test.y) / length(predicted.model2)))
## [1] "La precissió de l'arbre és del: 91.30 %"

Tenin una precissió de model força bona, de mes del 91% dels casos, però està una mica per sota que sense fer winnowed.

Fent servir la matriu de confusió veiem la qualitat del model també:

## creem la matriu per avaluar la prova
matrix.conf <- table(Class=predicted.model2,Predicted=test.y)
percent.correct <- 100 * sum(diag(matrix.conf)) / sum(matrix.conf)
print(sprintf("L'error de classificació és: %.2f %%",100 - percent.correct))
## [1] "L'error de classificació és: 8.70 %"
mosaicplot(matrix.conf)

Encara que la predicció sigui realment bona hi ha una cosa que no m’acaba d’agradar i és que en la predicció del ‘yes’ genera molt mes errors. Quan tenim un cas que ens diu que pot contractar un disposit bancari, té força errors. En canvi, el ‘no’, és molt precís. De totes maneres, és lleugerament millor que sense ‘winnow’.

confusionMatrix(predicted.model2,test.y, dnn = c("Prediction"))
## Confusion Matrix and Statistics
## 
##           NA
## Prediction    no   yes
##        no  10116   855
##        yes   143   362
##                                           
##                Accuracy : 0.913           
##                  95% CI : (0.9077, 0.9181)
##     No Information Rate : 0.894           
##     P-Value [Acc > NIR] : 4.875e-12       
##                                           
##                   Kappa : 0.382           
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.9861          
##             Specificity : 0.2975          
##          Pos Pred Value : 0.9221          
##          Neg Pred Value : 0.7168          
##              Prevalence : 0.8940          
##          Detection Rate : 0.8815          
##    Detection Prevalence : 0.9560          
##       Balanced Accuracy : 0.6418          
##                                           
##        'Positive' Class : no              
## 

Amb les estadístiques veiem que el ‘no’ té un 0.9221 i en canvi el ‘yes’ 0.7168 El valor de Kappa que tenim no és molt bo. Aquest valor és una altra forma d’interpretar l’aleatorietat explicat en el punt anterior. Ens ha baixat una mica fent servir menys atributs.

8.6 Estudi del model sense eliminar neteja de dades

Havíem reservat un model per a provar si millorava o empitjorava el fet d’eliminar outliers i alguns camps nulls.

bank.data <- read.csv('bank-discret.csv', sep = ',')

# separem el 70% de conjunt d'entrenament del 30% per a avaluar
random.data <- sample(1:nrow(bank.data), 0.70 * nrow(bank.data))
data.train <- bank.data[random.data,]
data.test <- bank.data[-random.data,]

# en les dades d'entrenament, separem les dades per a generar l'arbre de la dada objectiu (target)
train.x <- data.train[,1:16]
train.y <- data.train[,17]

# en les dades de test, separem les dades per a generar l'arbre de la dada objectiu (target)
test.x <- data.test[,1:16]
test.y <- data.test[,17]


model3 <- C50::C5.0(train.x, train.y, rules=TRUE, size=10,control = C5.0Control(winnow = TRUE))
predicted.model3 <- predict(model3, test.x, type="class")

matrix.conf <- table(Class=predicted.model3,Predicted=test.y)
mosaicplot(matrix.conf)

confusionMatrix(predicted.model3,test.y, dnn = c("Prediction"))
## Confusion Matrix and Statistics
## 
##           NA
## Prediction    no   yes
##        no  11716  1087
##        yes   263   498
##                                           
##                Accuracy : 0.9005          
##                  95% CI : (0.8953, 0.9055)
##     No Information Rate : 0.8831          
##     P-Value [Acc > NIR] : 7.382e-11       
##                                           
##                   Kappa : 0.3773          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.9780          
##             Specificity : 0.3142          
##          Pos Pred Value : 0.9151          
##          Neg Pred Value : 0.6544          
##              Prevalence : 0.8831          
##          Detection Rate : 0.8638          
##    Detection Prevalence : 0.9439          
##       Balanced Accuracy : 0.6461          
##                                           
##        'Positive' Class : no              
## 

La predicció en quant a accuracy és pitjor, el valor kappa (aletorietat) també és pitjor. La predicció del ‘yes’ ens cau a ‘0.6544’. Están ben netejades les dades perquè ens perjudica al classificar tenir valors nulls i outliers.

8.7 Reducció de la dimensionalitat (PCA)

Anem a reduïr encara més la dimensionalitat del dataset. Hi ha alguns atributs que potser ens fan embrutar l’arbre de decissió amb regles no necessaries

Primer veiem la matriu de correlació per veure quins atributs tenen impacte amb la variable de sortida target.

bank.data <- read.csv('bank-clean-discret.csv', sep = ',')
bank.data.tmp <- bank.data
for (colname in colnames(bank.data.tmp)) {
  bank.data.tmp[colname] <- lapply(bank.data.tmp[colname],as.integer)
}
corr.mat <- cor(bank.data.tmp)
# visualize it
library(corrplot)
## Warning: package 'corrplot' was built under R version 3.6.3
## corrplot 0.84 loaded
col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))
corrplot(corr.mat, method="color", col=col(200),   
         addCoef.col = "black", 
         tl.col="black", tl.srt=45,
         insig = "blank", 
         number.cex=0.5)

Veiem que els atributs que tenen més impacte en quant a la matriu de correlació serien (agafem > 0.07): education, balance, housing, contact, duration, pdays, previous

8.7.0.1 Estudi PCA (Principal Component Analysis)

En l’estudi de PCA el que es vol es trobar un nombre de dimensions menor de la mostra que expliquin el mateix que la mostra total. A nivell matemàtic funciona amb el concepte de ‘eigenvectors’ i ‘eigenvalues’, que són valors que simplifiquen una matriu dient el mateix que la matriu original: Eigenvector and Eigenvalue

El que fem és trobar la variànça total del conjunt de dades i intentarem reduir la dimensionalitat mantenint aquesta varinça al màxim. Per exemple, si tinguesim amb 10 dimensions un > 95% de la variança ja ens donariem per satisfets. L’elecció es realitza de manera que la primera component principal sigui la que major variància reculli; la segona ha de recollir la màxima variabilitat no recollida per la primera, i així successivament, triant un nombre que reculli un percentatge suficient de variància total.

Dividirem l’estudi en les dades quantitatives (PCA, Principal Component Analysis) i mixtes (FAMD, Factor Analysis of Mixed Data) per a dades qualitatives i quantitatives. Info

8.7.0.1.1 Dades quantitatives (PCA)
bank.pca <- read.csv('bank-clean-discret.csv', sep = ',')
quantitative <-  bank.pca[,c(1,6,10,12:15)]

# executem l'estudi PCA
pca <- prcomp(quantitative, scale. = TRUE)

# la funció mostra la variança de cada component 
fviz_eig(pca, choice = c("variance"), addlabels = TRUE)

# relació entre les dimensions i atributs.
fviz_contrib(pca, choice = "var", axes = 1:4)

# pve acumulada
PVE <- 100*pca$sdev^2/sum(pca$sdev^2)

plot(cumsum(PVE), type = "o", 
     ylab = "PVE acumulada", 
     xlab = "Componente principal", 
     col = "blue")

Totes les variables continues tenen una bona aportació encara que podriem reduïr el nombre ja que amb 4 dimensions estem aportant molta informació de la mostra.

En la tercera gràfica es pot veure que amb 5 dimensions tenim una variança acumulada de quasi un 80%. Podríem eliminar day, duration i campaign

## obetnim primer una mostra del total de forma aleatoria (hem guardat físicament aquest random per poder reproduïr el mateix i el carreguem de disc)
bank.data <- read.csv('bank-clean-discret.csv', sep = ',')
summary(bank.data)
##       age                 job           marital          education    
##  Min.   :1.000   blue-collar:8507   divorced: 4494   primary  : 5926  
##  1st Qu.:2.000   management :7896   married :22762   secondary:20981  
##  Median :3.000   technician :6615   single  :10997   tertiary :11346  
##  Mean   :3.046   admin.     :4565                                     
##  3rd Qu.:4.000   services   :3708                                     
##  Max.   :5.000   retired    :1492                                     
##                  (Other)    :5470                                     
##  default        balance        housing      loan            contact     
##  no :37485   Min.   :-1884.0   no :16220   no :31548   cellular :24919  
##  yes:  768   1st Qu.:   43.0   yes:22033   yes: 6705   telephone: 2169  
##              Median :  338.0                           unknown  :11165  
##              Mean   :  623.6                                            
##              3rd Qu.:  957.0                                            
##              Max.   : 3388.0                                            
##                                                                         
##       day            month          duration        campaign     
##  Min.   : 1.00   may    :12161   Min.   :1.000   Min.   : 1.000  
##  1st Qu.: 8.00   jul    : 6141   1st Qu.:2.000   1st Qu.: 1.000  
##  Median :16.00   aug    : 5320   Median :3.000   Median : 2.000  
##  Mean   :15.78   jun    : 4333   Mean   :3.003   Mean   : 2.776  
##  3rd Qu.:21.00   nov    : 2940   3rd Qu.:4.000   3rd Qu.: 3.000  
##  Max.   :31.00   apr    : 2450   Max.   :5.000   Max.   :58.000  
##                  (Other): 4908                                   
##      pdays          previous           poutcome       y        
##  Min.   : -1.0   Min.   :  0.0000   failure: 4099   no :34124  
##  1st Qu.: -1.0   1st Qu.:  0.0000   other  : 1551   yes: 4129  
##  Median : -1.0   Median :  0.0000   success: 1156              
##  Mean   : 40.3   Mean   :  0.5659   unknown:31447              
##  3rd Qu.: -1.0   3rd Qu.:  0.0000                              
##  Max.   :871.0   Max.   :275.0000                              
## 
bank.data.subset <- bank.data[,c("balance","age","previous","pdays","day","y")]
summary(bank.data.subset)
##     balance             age           previous            pdays      
##  Min.   :-1884.0   Min.   :1.000   Min.   :  0.0000   Min.   : -1.0  
##  1st Qu.:   43.0   1st Qu.:2.000   1st Qu.:  0.0000   1st Qu.: -1.0  
##  Median :  338.0   Median :3.000   Median :  0.0000   Median : -1.0  
##  Mean   :  623.6   Mean   :3.046   Mean   :  0.5659   Mean   : 40.3  
##  3rd Qu.:  957.0   3rd Qu.:4.000   3rd Qu.:  0.0000   3rd Qu.: -1.0  
##  Max.   : 3388.0   Max.   :5.000   Max.   :275.0000   Max.   :871.0  
##       day          y        
##  Min.   : 1.00   no :34124  
##  1st Qu.: 8.00   yes: 4129  
##  Median :16.00              
##  Mean   :15.78              
##  3rd Qu.:21.00              
##  Max.   :31.00
random.data <- read.csv('random-data.csv', sep = ',')
random.data <- unlist(random.data, use.names=FALSE)


data.train <- bank.data.subset[random.data,]
data.test <- bank.data.subset[-random.data,]

# en les dades d'entrenament, separem les dades per a generar l'arbre de la dada objectiu (target)
train.x <- data.train[,1:5]
train.y <- data.train[,6]

# en les dades de test, separem les dades per a generar l'arbre de la dada objectiu (target)
test.x <- data.test[,1:5]
test.y <- data.test[,6]

model31 <- C50::C5.0(train.x, train.y)
predicted.model31 <- predict(model31, test.x, type="class")


print(sprintf("La precissió de l'arbre és del: %.2f %%",100*sum(predicted.model31 == test.y) / length(predicted.model31)))
## [1] "La precissió de l'arbre és del: 89.85 %"
matrix.conf <- table(Class=predicted.model31,Predicted=test.y)
percent.correct <- 100 * sum(diag(matrix.conf)) / sum(matrix.conf)
print(sprintf("L'error de classificació és: %.2f %%",100 - percent.correct))
## [1] "L'error de classificació és: 10.15 %"
mosaicplot(matrix.conf)

confusionMatrix(predicted.model31,test.y, dnn = c("Prediction"))
## Confusion Matrix and Statistics
## 
##           NA
## Prediction    no   yes
##        no  10150  1056
##        yes   109   161
##                                          
##                Accuracy : 0.8985         
##                  95% CI : (0.8928, 0.904)
##     No Information Rate : 0.894          
##     P-Value [Acc > NIR] : 0.05853        
##                                          
##                   Kappa : 0.1852         
##                                          
##  Mcnemar's Test P-Value : < 2e-16        
##                                          
##             Sensitivity : 0.9894         
##             Specificity : 0.1323         
##          Pos Pred Value : 0.9058         
##          Neg Pred Value : 0.5963         
##              Prevalence : 0.8940         
##          Detection Rate : 0.8845         
##    Detection Prevalence : 0.9765         
##       Balanced Accuracy : 0.5608         
##                                          
##        'Positive' Class : no             
## 

La predicció no és gens bona i tenim un Kappa molt baix. Realment es troben a faltar les variables categòr

8.7.0.1.2 Dades mixtes (FAMD)

Com que el nostre dataset té dades categòriques i quantitatives, anem a fer l’estudi de dades mixtes amb la funció FAMD.

mixed <- read.csv('bank-clean-discret.csv', sep = ',')

res.famd <- FAMD(mixed, sup.var = 17, graph = FALSE, ncp = 40)

# relació entre les dimensions i atributs.
fviz_contrib(res.famd, choice = "var", axes = 1:40)

plot(res.famd$eig[,3], type = "o", 
     ylab = "PVE acumulada", 
     xlab = "Componente principal", 
     col = "blue")

Veiem en l’estudi que tenim 5 o 6 components que aporten moltíssim i la resta tenen un impacte similar. Anem a agafar un subconjunt del dataset que tenim guardat només amb els atributs: month, job, education, poutcome i contact.

## obetnim primer una mostra del total de forma aleatoria (hem guardat físicament aquest random per poder reproduïr el mateix i el carreguem de disc)
bank.data <- read.csv('bank-clean-discret.csv', sep = ',')
summary(bank.data)
##       age                 job           marital          education    
##  Min.   :1.000   blue-collar:8507   divorced: 4494   primary  : 5926  
##  1st Qu.:2.000   management :7896   married :22762   secondary:20981  
##  Median :3.000   technician :6615   single  :10997   tertiary :11346  
##  Mean   :3.046   admin.     :4565                                     
##  3rd Qu.:4.000   services   :3708                                     
##  Max.   :5.000   retired    :1492                                     
##                  (Other)    :5470                                     
##  default        balance        housing      loan            contact     
##  no :37485   Min.   :-1884.0   no :16220   no :31548   cellular :24919  
##  yes:  768   1st Qu.:   43.0   yes:22033   yes: 6705   telephone: 2169  
##              Median :  338.0                           unknown  :11165  
##              Mean   :  623.6                                            
##              3rd Qu.:  957.0                                            
##              Max.   : 3388.0                                            
##                                                                         
##       day            month          duration        campaign     
##  Min.   : 1.00   may    :12161   Min.   :1.000   Min.   : 1.000  
##  1st Qu.: 8.00   jul    : 6141   1st Qu.:2.000   1st Qu.: 1.000  
##  Median :16.00   aug    : 5320   Median :3.000   Median : 2.000  
##  Mean   :15.78   jun    : 4333   Mean   :3.003   Mean   : 2.776  
##  3rd Qu.:21.00   nov    : 2940   3rd Qu.:4.000   3rd Qu.: 3.000  
##  Max.   :31.00   apr    : 2450   Max.   :5.000   Max.   :58.000  
##                  (Other): 4908                                   
##      pdays          previous           poutcome       y        
##  Min.   : -1.0   Min.   :  0.0000   failure: 4099   no :34124  
##  1st Qu.: -1.0   1st Qu.:  0.0000   other  : 1551   yes: 4129  
##  Median : -1.0   Median :  0.0000   success: 1156              
##  Mean   : 40.3   Mean   :  0.5659   unknown:31447              
##  3rd Qu.: -1.0   3rd Qu.:  0.0000                              
##  Max.   :871.0   Max.   :275.0000                              
## 
bank.data.subset <- bank.data[,c("month","job","education","poutcome","contact","y")]
summary(bank.data.subset)
##      month                job           education        poutcome    
##  may    :12161   blue-collar:8507   primary  : 5926   failure: 4099  
##  jul    : 6141   management :7896   secondary:20981   other  : 1551  
##  aug    : 5320   technician :6615   tertiary :11346   success: 1156  
##  jun    : 4333   admin.     :4565                     unknown:31447  
##  nov    : 2940   services   :3708                                    
##  apr    : 2450   retired    :1492                                    
##  (Other): 4908   (Other)    :5470                                    
##       contact        y        
##  cellular :24919   no :34124  
##  telephone: 2169   yes: 4129  
##  unknown  :11165              
##                               
##                               
##                               
## 
random.data <- read.csv('random-data.csv', sep = ',')
random.data <- unlist(random.data, use.names=FALSE)


data.train <- bank.data.subset[random.data,]
data.test <- bank.data.subset[-random.data,]

# en les dades d'entrenament, separem les dades per a generar l'arbre de la dada objectiu (target)
train.x <- data.train[,1:5]
train.y <- data.train[,6]

# en les dades de test, separem les dades per a generar l'arbre de la dada objectiu (target)
test.x <- data.test[,1:5]
test.y <- data.test[,6]

model4 <- C50::C5.0(train.x, train.y, rules=TRUE)
summary(model4)
## 
## Call:
## C5.0.default(x = train.x, y = train.y, rules = TRUE)
## 
## 
## C5.0 [Release 2.07 GPL Edition]      Tue Jun 16 22:30:50 2020
## -------------------------------
## 
## Class specified by attribute `outcome'
## 
## Read 26777 cases (6 attributes) from undefined.data
## 
## Rules:
## 
## Rule 1: (25949/2393, lift 1.0)
##  poutcome in {failure, other, unknown}
##  ->  class no  [0.908]
## 
## Rule 2: (828/309, lift 5.8)
##  poutcome = success
##  ->  class yes  [0.627]
## 
## Default class: no
## 
## 
## Evaluation on training data (26777 cases):
## 
##          Rules     
##    ----------------
##      No      Errors
## 
##       2 2702(10.1%)   <<
## 
## 
##     (a)   (b)    <-classified as
##    ----  ----
##   23556   309    (a): class no
##    2393   519    (b): class yes
## 
## 
##  Attribute usage:
## 
##  100.00% poutcome
## 
## 
## Time: 0.1 secs
predicted.model4 <- predict(model4, test.x, type="class")
print(sprintf("La precissió de l'arbre és del: %.2f %%",100*sum(predicted.model4 == test.y) / length(predicted.model4)))
## [1] "La precissió de l'arbre és del: 90.39 %"

Obtenim un arbre molt més simple amb unes regles també més simples però el valor de la predicció ens ha baixar una mica.

predicted.model4 <- predict(model4, test.x, type="class")
print(sprintf("La precissió de l'arbre és del: %.2f %%",100*sum(predicted.model4 == test.y) / length(predicted.model4)))
## [1] "La precissió de l'arbre és del: 90.39 %"

Amb una predicció lleugerament inferior al model inicial, però no moltíssim.

matrix.conf <- table(Class=predicted.model4,Predicted=test.y)
percent.correct <- 100 * sum(diag(matrix.conf)) / sum(matrix.conf)
print(sprintf("L'error de classificació és: %.2f %%",100 - percent.correct))
## [1] "L'error de classificació és: 9.61 %"
mosaicplot(matrix.conf)

confusionMatrix(predicted.model4,test.y, dnn = c("Prediction"))
## Confusion Matrix and Statistics
## 
##           NA
## Prediction    no   yes
##        no  10152   996
##        yes   107   221
##                                           
##                Accuracy : 0.9039          
##                  95% CI : (0.8983, 0.9092)
##     No Information Rate : 0.894           
##     P-Value [Acc > NIR] : 0.000245        
##                                           
##                   Kappa : 0.2524          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.9896          
##             Specificity : 0.1816          
##          Pos Pred Value : 0.9107          
##          Neg Pred Value : 0.6738          
##              Prevalence : 0.8940          
##          Detection Rate : 0.8846          
##    Detection Prevalence : 0.9714          
##       Balanced Accuracy : 0.5856          
##                                           
##        'Positive' Class : no              
## 

El que si que ens ha perjudicat és en el valor de Kappa, on sembla que l’aletorietat és bastant pitjor

model4 <- C50::C5.0(train.x, train.y, rules=FALSE)
plot(model4)


9 Classificació: Hi ha hagut millora en capacitat predictiva, després d’aplicar PCA/SVD? A què creus que és degut?


Hem fet 4 estudis:

Conclusions: - Podem dir que la classificació amb PCA(winnowed) i sense PCA ens ha donat valors realment molt similars. - Amb PCA (winnowed), hem millorat en la predicció general a causa que tenim millor predicció del ‘yes’. - Sense PCA hem obtingut bons valors i el millor valor de Kappa. - La millora en la predicció aplicant PCA(winnowed) és molt petita. Això pot venir donat per el fet que necessitem moltes dimensions per a tenir una variança acumulada per sobre de 95%. No podem desfer-nos de moltes varibles.


10 Classificació: Altres tipus d’algoritmes


10.1 Random Forest

Un Random Forest és un conjunt d’arbres de decisió combinats amb bagging. A l’usar bagging, el que en realitat està passant, és que diferents arbres veuen diferents porcions de les dades. Cap arbre veu totes les dades d’entrenament. Això fa que cada arbre s’entreni amb diferents mostres de dades per a un mateix problema. D’aquesta manera, a l’combinar els seus resultats, uns errors es compensen amb altres i tenim una predicció que generalitza millor.

bank.data <- read.csv('bank-clean-discret.csv', sep = ',')

random.data <- read.csv('random-data.csv', sep = ',')
random.data <- unlist(random.data, use.names=FALSE)


train <- bank.data[random.data,]
test <- bank.data[-random.data,]

test.x <- test[,1:16]
test.y <- test[,17]

train_rf <- randomForest(y~.,data = train, ntree = 500)

predicted.model5 <- predict(train_rf,test.x)


print(sprintf("La precissió de l'arbre és del: %.2f %%",100*sum(predicted.model5 == test.y) / length(predicted.model5)))
## [1] "La precissió de l'arbre és del: 91.08 %"
matrix.conf <- table(Class=predicted.model5,Predicted=test.y)
percent.correct <- 100 * sum(diag(matrix.conf)) / sum(matrix.conf)
print(sprintf("L'error de classificació és: %.2f %%",100 - percent.correct))
## [1] "L'error de classificació és: 8.92 %"
mosaicplot(matrix.conf)

confusionMatrix(predicted.model5,test.y, dnn = c("Prediction"))
## Confusion Matrix and Statistics
## 
##           NA
## Prediction    no   yes
##        no  10033   798
##        yes   226   419
##                                           
##                Accuracy : 0.9108          
##                  95% CI : (0.9054, 0.9159)
##     No Information Rate : 0.894           
##     P-Value [Acc > NIR] : 1.13e-09        
##                                           
##                   Kappa : 0.4064          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.9780          
##             Specificity : 0.3443          
##          Pos Pred Value : 0.9263          
##          Neg Pred Value : 0.6496          
##              Prevalence : 0.8940          
##          Detection Rate : 0.8743          
##    Detection Prevalence : 0.9438          
##       Balanced Accuracy : 0.6611          
##                                           
##        'Positive' Class : no              
## 
  • Random Forest:

    • Accuracy: 0.9108
    • Kappa: 0.4064
    • Prediction ‘no’: 0.9263
    • Prediction ‘yes’: 0.6496

Tenim bona predicció i Kappa però la predicció del ‘no’ és força baixa (0.6496)

10.2 CART

A diferència del ID3, aquest tipus d’arbre:

  • Mesura d’homogeneïtat i criteri d’aturada en el procés de partició i divisió de l’arbre a partir de l’índex de Gini, encara que hi ha variants que n’escullen d’altres. Aquesta mesura estableix el millor separador com el que redueix la diversitat de les diferents particions obtingudes (subarbres).

  • La poda es realitza de la següent forma:

    • Generar diversos subarbres podats “interessants”.
    • Obtenir estimacions de l’error de cadascun d’aquests subarbres.
    • Escollir el subarbre que presenti la millor estimació.
  • Són binaris; a cada node hi ha un punt de tall (per un procediment semblant al que s’ha explicat per a trobar punts de tall en la discretització d’atributs continus) que separa en dos el conjunt d’observacions.

  • L’algorisme CART pot treballar amb atributs continus (tot i que les modificacions de ID3 també ho poden fer com hem pogut veure en l’exemple anterior C4.5).

  • Pot fer tant classificació com regressió: en el primer cas, la variable per a predir ha de ser categòrica amb un valor per a cada classe possible.

bank.data <- read.csv('bank-clean-discret.csv', sep = ',')
random.data <- read.csv('random-data.csv', sep = ',')
random.data <- unlist(random.data, use.names=FALSE)

train <- bank.data[random.data,]
test <- bank.data[-random.data,]

test.x <- test[,1:16]
test.y <- test[,17]

heart.tree <- rpart(y ~ .,method="class", data=train)
summary(heart.tree)
## Call:
## rpart(formula = y ~ ., data = train, method = "class")
##   n= 26777 
## 
##           CP nsplit rel error    xerror       xstd
## 1 0.02403846      0 1.0000000 1.0000000 0.01749460
## 2 0.01236264      4 0.9007555 0.9007555 0.01670403
## 3 0.01000000      6 0.8760302 0.8839286 0.01656405
## 
## Variable importance
## duration poutcome    month  contact 
##       53       36        6        5 
## 
## Node number 1: 26777 observations,    complexity param=0.02403846
##   predicted class=no   expected loss=0.10875  P(node) =1
##     class counts: 23865  2912
##    probabilities: 0.891 0.109 
##   left son=2 (21471 obs) right son=3 (5306 obs)
##   Primary splits:
##       duration < 4.5    to the left,    improve=597.0337, (0 missing)
##       poutcome splits as  LLRL,         improve=458.6320, (0 missing)
##       month    splits as  LLRLLLLRLLRR, improve=270.8700, (0 missing)
##       pdays    < 8.5    to the left,    improve=131.0225, (0 missing)
##       previous < 0.5    to the left,    improve=129.9314, (0 missing)
##   Surrogate splits:
##       balance  < 3386.5 to the left,  agree=0.802, adj=0, (0 split)
##       previous < 53     to the left,  agree=0.802, adj=0, (0 split)
## 
## Node number 2: 21471 observations,    complexity param=0.02403846
##   predicted class=no   expected loss=0.05626193  P(node) =0.8018449
##     class counts: 20263  1208
##    probabilities: 0.944 0.056 
##   left son=4 (20854 obs) right son=5 (617 obs)
##   Primary splits:
##       poutcome splits as  LLRL,         improve=342.36130, (0 missing)
##       month    splits as  LLRLLLLRLLRR, improve=216.18880, (0 missing)
##       pdays    < 0      to the left,    improve= 98.61017, (0 missing)
##       previous < 0.5    to the left,    improve= 98.61017, (0 missing)
##       housing  splits as  RL,           improve= 53.97091, (0 missing)
## 
## Node number 3: 5306 observations,    complexity param=0.02403846
##   predicted class=no   expected loss=0.3211459  P(node) =0.1981551
##     class counts:  3602  1704
##    probabilities: 0.679 0.321 
##   left son=6 (5095 obs) right son=7 (211 obs)
##   Primary splits:
##       poutcome splits as  LLRL,         improve=91.42517, (0 missing)
##       contact  splits as  RRL,          improve=72.66632, (0 missing)
##       month    splits as  LLRLLLLRLLRR, improve=48.09435, (0 missing)
##       housing  splits as  RL,           improve=39.70613, (0 missing)
##       pdays    < 28.5   to the left,    improve=36.09150, (0 missing)
## 
## Node number 4: 20854 observations
##   predicted class=no   expected loss=0.04090342  P(node) =0.7788027
##     class counts: 20001   853
##    probabilities: 0.959 0.041 
## 
## Node number 5: 617 observations,    complexity param=0.02403846
##   predicted class=yes  expected loss=0.4246353  P(node) =0.02304216
##     class counts:   262   355
##    probabilities: 0.425 0.575 
##   left son=10 (177 obs) right son=11 (440 obs)
##   Primary splits:
##       duration < 2.5    to the left,    improve=44.239210, (0 missing)
##       month    splits as  RRRRLRRRLRRR, improve=13.032810, (0 missing)
##       housing  splits as  RL,           improve=10.480080, (0 missing)
##       job      splits as  LLLRRRRRRLR,  improve= 5.212144, (0 missing)
##       pdays    < 199.5  to the right,   improve= 5.145055, (0 missing)
##   Surrogate splits:
##       contact  splits as  RRL,        agree=0.724, adj=0.040, (0 split)
##       campaign < 6.5    to the right, agree=0.720, adj=0.023, (0 split)
##       pdays    < 606    to the right, agree=0.716, adj=0.011, (0 split)
##       default  splits as  RL,         agree=0.715, adj=0.006, (0 split)
## 
## Node number 6: 5095 observations,    complexity param=0.01236264
##   predicted class=no   expected loss=0.3022571  P(node) =0.1902752
##     class counts:  3555  1540
##    probabilities: 0.698 0.302 
##   left son=12 (1500 obs) right son=13 (3595 obs)
##   Primary splits:
##       contact splits as  RRL,          improve=54.21724, (0 missing)
##       month   splits as  LLRLLLLRLLRR, improve=36.53465, (0 missing)
##       housing splits as  RL,           improve=28.85552, (0 missing)
##       pdays   < 371.5  to the left,    improve=16.16480, (0 missing)
##       job     splits as  RLLLRRRLRRR,  improve=14.23330, (0 missing)
##   Surrogate splits:
##       month    splits as  RRRRRRLRLRRR, agree=0.848, adj=0.483, (0 split)
##       campaign < 24.5   to the right,   agree=0.706, adj=0.002, (0 split)
##       balance  < -891.5 to the left,    agree=0.706, adj=0.001, (0 split)
## 
## Node number 7: 211 observations
##   predicted class=yes  expected loss=0.2227488  P(node) =0.007879897
##     class counts:    47   164
##    probabilities: 0.223 0.777 
## 
## Node number 10: 177 observations
##   predicted class=no   expected loss=0.2768362  P(node) =0.006610151
##     class counts:   128    49
##    probabilities: 0.723 0.277 
## 
## Node number 11: 440 observations
##   predicted class=yes  expected loss=0.3045455  P(node) =0.01643201
##     class counts:   134   306
##    probabilities: 0.305 0.695 
## 
## Node number 12: 1500 observations
##   predicted class=no   expected loss=0.1893333  P(node) =0.05601822
##     class counts:  1216   284
##    probabilities: 0.811 0.189 
## 
## Node number 13: 3595 observations,    complexity param=0.01236264
##   predicted class=no   expected loss=0.3493741  P(node) =0.134257
##     class counts:  2339  1256
##    probabilities: 0.651 0.349 
##   left son=26 (3347 obs) right son=27 (248 obs)
##   Primary splits:
##       month   splits as  LLRLLLRRLLRR, improve=46.61047, (0 missing)
##       housing splits as  RL,           improve=18.72056, (0 missing)
##       pdays   < 371.5  to the left,    improve=12.98456, (0 missing)
##       job     splits as  RLLLRRRRRRR,  improve=10.95079, (0 missing)
##       balance < 1597.5 to the left,    improve=10.86097, (0 missing)
##   Surrogate splits:
##       day   < 1.5    to the right, agree=0.933, adj=0.028, (0 split)
##       pdays < 472.5  to the left,  agree=0.932, adj=0.008, (0 split)
## 
## Node number 26: 3347 observations
##   predicted class=no   expected loss=0.3274574  P(node) =0.1249953
##     class counts:  2251  1096
##    probabilities: 0.673 0.327 
## 
## Node number 27: 248 observations
##   predicted class=yes  expected loss=0.3548387  P(node) =0.00926168
##     class counts:    88   160
##    probabilities: 0.355 0.645
predicted.model6 <- predict(heart.tree, newdata = test.x, type = "class")
matrix.conf <- table(Class=test.y,Predicted=predicted.model6)
percent.correct <- 100 * sum(diag(matrix.conf)) / sum(matrix.conf)
print(sprintf("L'error de classificació és: %.2f %%",100 - percent.correct))
## [1] "L'error de classificació és: 9.15 %"
mosaicplot(matrix.conf)

confusionMatrix(predicted.model6,test.y, dnn = c("Prediction"))
## Confusion Matrix and Statistics
## 
##           NA
## Prediction    no   yes
##        no  10153   944
##        yes   106   273
##                                           
##                Accuracy : 0.9085          
##                  95% CI : (0.9031, 0.9137)
##     No Information Rate : 0.894           
##     P-Value [Acc > NIR] : 1.29e-07        
##                                           
##                   Kappa : 0.3072          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.9897          
##             Specificity : 0.2243          
##          Pos Pred Value : 0.9149          
##          Neg Pred Value : 0.7203          
##              Prevalence : 0.8940          
##          Detection Rate : 0.8847          
##    Detection Prevalence : 0.9670          
##       Balanced Accuracy : 0.6070          
##                                           
##        'Positive' Class : no              
## 

Obtenim una predicció molt similar, quasi del 0.91 amb un Kappa del 0.3072. La predicció del ‘no’ és molt bona, 0.9149 i em millorat la del yes a un valor de 0.7203

  • Cart:

    • Accuracy: 0.91
    • Kappa: 0.3072.
    • Prediction ‘no’: 0.9149
    • Prediction ‘yes’: 0.7203
fancyRpartPlot(heart.tree, caption = NULL)

Aquest arbre és molt més interpretable que el de l’algoritme C5.0. CART ha realitzat l’eliminació d’atributs i només ens hem quedat amb duration, poutcome, month i contact.

  • És molt important que les trucades tinguin una durada superior a 4.5 minuts.
  • El fet que una campanya anterior hagués estat success o no és el següent camp en importància.

10.3 k-Nearest Neighbors

Per finalitzar farem servir un altre mètode per classificar basat en la teoria que vam veure a clustering. El que fem és trobar valors veïns similars a partir de les distàncies dels diferents atributs, és a dir, troba similaritat de posició amb N dimensions. Hem de definir un valor de K. Aquest valor indica el nombre de veïns per a definir un nou grup o en aquest cas predir el resultat.

  • Farem servir la funció de RStudio (knn).
  • També normalitzarem els atributs per calcular les distancies.
  • Agafarem només els valors numèrics. Podem agafar els categòrics i passar-los a numèrics també.

Preparem les dades:

# Knn
bank.data <- read.csv('bank-clean-discret.csv', sep = ',')

bank.data <- bank.data[,c("poutcome","month","pdays","day","housing","duration","contact","loan","balance","y")]

# separem el 70% de conjunt d'entrenament del 30% per a avaluar
random.data <- read.csv('random-data.csv', sep = ',')
random.data <- unlist(random.data, use.names=FALSE)

train <- bank.data[random.data,]
test <- bank.data[-random.data,]

# en les dades d'entrenament, separem les dades per a generar l'arbre de la dada objectiu (target)
train.x <- train[,1:9]
train.y <- train[,10]

# en les dades de test, separem les dades per a generar l'arbre de la dada objectiu (target)
test.x <- test[,1:9]
test.y <- test[,10]

## creem la funció per normalitzar
ni <-function (x) {(x -min (x)) / (max (x) -min (x))} 

# passem a numèric i normalitzem tant test(x,y) com train(x,y)
train.x.norm <- train.x
train.x.norm <- lapply(train.x.norm[1:9],as.integer)
train.x.norm <- as.data.frame(lapply(train.x.norm[1:9],ni))

train.y.norm <- train.y
train.y.norm <- as.numeric(train.y.norm)
train.y.norm[train.y.norm == 1] <- 0
train.y.norm[train.y.norm == 2] <- 1

test.x.norm <- test.x
test.x.norm[1:9] <- lapply(test.x.norm[1:9],as.integer)
test.x.norm[1:9] <- as.data.frame(lapply(test.x.norm[1:9],ni))

test.y.norm <- test.y
test.y.norm <- as.numeric(test.y.norm)
test.y.norm[test.y.norm == 1] <- 0
test.y.norm[test.y.norm == 2] <- 1

Calculem un valor optim de k. Veins propers.

 i=1                          # declaration to initiate for loop
 k.optm=1                     # declaration to initiate for loop
 for (i in 1:75){ 
     predicted.model7 <-  knn(train=train.x.norm, test=test.x.norm, cl=train.y.norm, k=i)
     k.optm[i] <- 100*sum(predicted.model7 == test.y.norm) / length(predicted.model7)
     k=i
 }
 plot(k.optm, type="b", xlab="K- Value",ylab="Accuracy level") 

Veint la gràfica optem per a fer servir un k entre 20 i 40

predicted.model7 <- knn(train.x.norm,test.x.norm,cl=train.y.norm,k=25)
print(sprintf("La precissió de l'arbre és del: %.2f %%",100*sum(predicted.model7 == test.y.norm) / length(predicted.model7)))
## [1] "La precissió de l'arbre és del: 90.59 %"
matrix.conf <- table(Class=predicted.model7,Predicted=test.y)
percent.correct <- 100 * sum(diag(matrix.conf)) / sum(matrix.conf)
print(sprintf("L'error de classificació és: %.2f %%",100 - percent.correct))
## [1] "L'error de classificació és: 9.41 %"
mosaicplot(matrix.conf)

confusionMatrix(predicted.model7,as.factor(test.y.norm), dnn = c("Prediction"))
## Confusion Matrix and Statistics
## 
##           NA
## Prediction     0     1
##          0 10108   929
##          1   151   288
##                                           
##                Accuracy : 0.9059          
##                  95% CI : (0.9004, 0.9112)
##     No Information Rate : 0.894           
##     P-Value [Acc > NIR] : 1.3e-05         
##                                           
##                   Kappa : 0.309           
##                                           
##  Mcnemar's Test P-Value : < 2e-16         
##                                           
##             Sensitivity : 0.9853          
##             Specificity : 0.2366          
##          Pos Pred Value : 0.9158          
##          Neg Pred Value : 0.6560          
##              Prevalence : 0.8940          
##          Detection Rate : 0.8808          
##    Detection Prevalence : 0.9617          
##       Balanced Accuracy : 0.6110          
##                                           
##        'Positive' Class : 0               
## 
  • k-Nearest Neighbors:

    • Accuracy: 0.9059
    • Kappa: 0.309
    • Prediction ‘no’: 0.9158
    • Prediction ‘yes’: 0.6560

11 Classificació: Més Conclusions